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INTRODUCTION 

In  recent  years,  engineers  have  become  more  and  more  concerned  with 
optimization.     One  of  the  reasons  for  this  concern  is  that  intensive  compe- 
tition in  the  industries  in  general  and  in  the  chemical  process  industries  in 
particular  makes   it  more  necessary  than  ever  that  equipment  and  systems  be 
designed  and  operated  at  peak  performance.     Even  marginal  savings  can  be 
extremely  vital  in  this  competitive  environment  (1). 

A  simple  optimization  problem  corresponds  to  the  findings  of  the  extreme 
value  of  a  function  in  calculus   (2).     This  can  be  accomplished  by  differentia- 
tion.    As  the  problem  becomes  more  involved,  partial  differentiation  and  the 
calculus  of  variation  (3)  may  have  to  be  used.     More  often  than  not,   optimal 
problems  in  engineering  and  industry  cannot  be  solved  by  direct  applications 
of  these  conventional  mathematical  methods.     A  variety  of  approaches  more 
sophisticated  than  the  conventional  methods  has  been  proposed  to  solve  com- 
plex problems.     Among  them  are  Dynamic  Programming  (h)  and  the  Maximum  Prin- 
ciple  (5). 

Using  the  method  of  dynamic  programming,  Aris  has  solved  a  number  of 
problems  dealing  with  chemical  reactor  design  (6);  Aris  and  co-workers  have 
solved  a  cross-current  extraction  problem  (7);   Kudd  has  investigated  a 
reliability  problem  in  chemical  system  design  (8)  and  the  optimal  use  of 
limited  resources   (9).     Other  applications  of  dynamic  programming  to  chemical 
engineering  problems  can  be  found  in  the  works  of  Roberts   (10,11),  Dranoff 
et  al.   (12),  and  Mitten  et  al.   (13).     Katz's  algorithm  of  the  discrete  maxi- 
mum principle  has  been  successfully  applied  to  several  chemical  processes  in 
the  works  of  Ahn  et  al.    (lU)  and  Wang  (15). 


In  this  report  a  brief  review  of  these  two  approaches  and  their  applica- 
tions to  several  chemical  reactor  systems  is  presented. 

DYNAMIC  PROGRAMMING 

The  basic  notion  of  the  principle  of  optimality  from  which  dynamic  pro- 
gramming algorithm  is  derived  is  introduced  by  considering  a  simple  example. 
The  problem  is  formulated  according  to  the  dynamic  programming  algorithm.  Two 
techniques  frequently  used  to  reduce  complexities  of  problems  —  the  use  of 
the  Lagrange  Multiplier  and  the  ratio  of  the  state  variables  —  are  then  in- 
troduced by  again  using  examples. 

Basic  Problem — Minimum  Length  of  Directed  Network 

In  Figure  1,  the  circles  represent  the  nodes,  and  the  numbers  on  the 
lines  which  connect  nodes  denote  the  distances  between  the  nodes.  Suppose  that 
the  problem  is  to  find  the  shortest  path  from  node  1  to  node  6  with  an  addi- 
tional restriction  that  the  passages  of  a  particle  should  be  allowed  along  any 
path  only  in  the  directions  of  increasing  nodal  index  (16).  Common  sense  tells 
us  that  we  can  enumerate  all  possible  paths  between  node  1  and  node  6,  and  then 
pick  the  shortest  one.  For  convenience  the  number  in  the  parentheses  is  used 
to  denote  the  node  number. 

path  (D—^—OO— (6)  -  13 

path  (1)_(2)_ (S)_(6)  =  7 

path  (1)_(2)_(U  )—{$)— (6)  -  13 

path  (1  )_(2 )— (3  )_($)— (6)  =  8 

path  (1)_(2)_(3)_(10— (6)  -  10 

path  (1)_(2)_(3)_(U)_ (5)_(6)  =  10 
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path  (l)_(3)_(li)_(6)  -  10 

path  (1)_(3)_00— (5)— (6)  -  10 

path  (1)_(3)_(5)_(6)  -  8 

The  shortest  path  is  the  second  one  listed,  i.e.,  (l) (2) (£)_(6)  whose 

total  length  is  7. 

But  enumeration  of  all  the  possibilities  would  be  prohibitive  in  a 
problem  of  large  size,  employing  dynamic  programming  treatment,  we  work 
backwards  as  follows: 

from  (5)  to  (6)  =  2 

from  (It)  to  (6)  -  5  either  by  way  of  (k)   to  (6)  directly  or  (h)—($)—(6). 
To  go  from  (3)  to  (6)  there  are  two  choices,  either  by  way  of  (k)   or  (5).  If 
we  proceed  by  way  of  (h),   the  distance  from  (3)  to  (6)  is  the  sum  of  (3)  to 
(It)  and  (h)   to  (6),  that  is,  2  +  5  =  7.  This  is  true  for  both  paths: 
(3)— (I;)— (6)  or  (3 )—(!*)_ (5)— (6).  If  we  proceed  by  way  of  (5),  we  have: 
(3)— (5)— (6)  -  5.  Thus,  the  shortest  path  from  (3)  to  (6)  is: 
(3)—  (5)— (6)  »  5.  By  the  same  argument  there  are  three  alternatives  lead- 
ing from  (2)  to  (6),  namely,  (2) — (U)_(6)  =  12,  (2)_(3)_(S)_(6)  =  7 
and  (2)_(5)_(6)  -  6.  It  follows  that  the  shortest  path  from  (2)  to  (6) 
(2)_(5)_(6)  =  6.  Starting  from  (1),  we  can  go  either  by  way  of  (2)  or  (3) 
to  reach  (6).  The  distance  from  (1)  to  (3)  is  3  and  the  shortest  distance  from 
(3)  to  (6)  is  (3)_(5)_ (6)  -  J.  As  soon  as  the  particle  reaches  (3),  it 
must  follow  the  shortest  path  from  (3)  leading  to  (6)  in  order  to  make  the 
overall  path  the  shortest.  Therefore,  the  shortest  path  from  (1)  to  (6)  by 
»ay  of  (3)  is:  (1)_(3)_(5)_(6)  .  8  and  from  (1)  to  (6)  by  way  of  (2)  is 
(1)_(2)_(5)_(6)  -  7.  The  choice  of  the  shortest  path  from  (1)  to  (6) 
will  naturally  be  (1)_(2)_(5)_ (6)  -  7. 


The  main  difference  in  the  two  methods  is  that  in  the  first,  the  process 
of  enumerating  all  the  paths,  we  do  not  have  any  idea  which  path  will  eventu- 
ally be  the  answer.  We  are  only  exhausting  all  the  possible  cases.  However, 
the  second  method  allows  us  to  leave  out  all  the  impossible  cases.  For  example, 
(2) — (S) — (6)  -  6  is  always  the  shortest  path  between  (2)  and  (6)  and,  there- 
fore, every  time  we  proceed  by  way  of  (2),  we  discard  all  other  paths  branch- 
ing out  from  (2). 

What  is  illustrated  above  represents  a  very  simple  version  of  Bellman's 
Principle  of  Optimality  (ii)  which  states  that  "An  optimal  policy  has  the 
property  that  whatever  the  initial  state  and  initial  decision::,  are  the  remain- 
ing decisions  must  constitute  an  optimal  polity  with  regard  to  the  state  re- 
sulting from  the  first  decision."     Dynamic  programming  is  the  algorithm  based 
on  this  principle.     The  advantage  of  this  approach  is  to  reduce  a  problem 
with  large  dimensions  to  a  problem  of  solving  simultaneous  equations  each 
with  small  dimensions.     Then  the  use  of  a  high  speed  computer  becomes  partic- 
ularly suitable. 

Before  describing  in  detail  the  formal  algorithm  of  dynamic  programming, 
a  basic  terminology  should  be  presented  (9,1$). 

1.  State,  x:     The  condition  of  the  entry  to  be  transformed  in  the  process; 
for  instance,   the  concentration  of  the  reacting  substances  in  the 
reactor,  or  the  position  of  the  node  at  which  the  particle  is  currently 
situated. 

2.  Decision,  or  control  variable,   9:     A  decision  or  control  variable  is 
defined  as  a  way  of  manipulation  made  to  transform  the  entries  from  one 
state  to  another;  for  exanple,   the  temperature  or  pressure  in  a  reactor, 
or  the  decision  of  the  path  to  follow  in  the  net-work  problem. 


3.   Transformation,  T:  The  transformation  operator  describes  the  way  in 

which  a  given  decision  or  control  action  transforms  the  entries  from  one 
state  to  another,  for  example;  the  kinetic  equations  or  material  balance 
equations  and  others. 

h.       Interval  profit,  P:  The  interval  profit  is  the  profit  or  any  physical 
quantity  resulting  from  the  transformation  of  the  entries  in  the  interval 
time  or  space  between  two  decisions.  It  is  a  function  of  the  state  and 
the  decision  made.  An  example,  is  the  increase  in  concentration  of  a 
desired  product  in  one  backmix  chemical  reactor  of  a  multistage  process. 
Another  example  is  the  distance  between  two  nodes  in  the  interval. 

5.  Stage,  n:  This  represents  the  typical  unit  of  a  process  or  the  typical 
time  periods,  or  any  real  or  abstract  notion  of  intervals  in  which  cer- 
transformations  take  place. 

6.  Objective  function,  0:  This  is  summation  of  the  interval  profits,  it  is 
the  net  profit  of  the  entire  process j  for  example,  the  total  distance 

of  the  paths  from  node  1  to  6.  Another  example  is  the  change  in  the 
concentration  of  the  desired  product  between  the  first  and  last  stages 
of  a  multistage  process. 

The  optimal  sequence  of  decisions  is  defined  as  that  sequence  of 
decisions  which  maximizes  the  objective  function. 

7.  Constraints:  These  are  the  limitations  imposed  on  the  decisions; 
examples  are,  only  the  path  in  the  direction  of  increasing  nodal  index 
is  permitted  in  the  network  problem  or  the  limitations  on  the  tempera- 
ture that  may  be  imposed  in  designing  a  reactor,  or  the  limitations  of 
relative  cost  of  the  product. 


According  to  the  dynamic  programming  algorithm  the  net-work  problem  may- 
be reformulated  as  follows   (16): 

Let  £j_  be  the  shortest  path  from  node  i  to  6.     Pj»  is  the  distance  from 
node  i  to  j,  and  the  constraints  demand  that  j>i  and  ?^  ""^when  node  i  and 
0  are  not  connected.     Also,   it  can  be  noted  that  the  distance  from  node  6 
to  itself  is  zero. 

°6  "  0  (1) 

In  general 

ft-g(Pjj*fd)  -  (2) 

To  find  the  optimal  objective  function,   in  this  case  the  shortest  distance, 
+ 
0*  ■  t.f  we  employ  the  equation   (2)  repeatedly, 

fj-Mta     (PSj+fj).Min(P56+f6).2  +  0-2 

This  is  to  say  that,  if  the  system  were  in  state  5,  an  optimal  continuation 
would  be  to  proceed  to  node  6,  with  a  corresponding  increase  in  length  of 
2  units.  This  continuation  is  designated  by  Jc  ■  6. 

This  shows  that  both  paths  are  equal  in  distance  and,  therefore,  J  =  5  or  6. 

h 


The  star,  *,  indicates  optimum  conditions  in  general.  For  instance,  8*  is 
the  optimal  decision  and  0*  is  the  optimum  profit,  i.e.  the  objective  func- 
tion. 


/p  +  f 


8-Jj2;ji);«-(i:!)  =6  Vs 


25   5, 


,1  +  6, 


'1   2 


0  -f  =  Min  (    +f  )  -Min  C,  .")  -  7,  J 

Thus  the  objective  function  or  profit  is  7.  '  The  track  of  path  will  be 
i±   =  2,  j2  =  5,  j  ■  6,  i.e.,  (1) — (2) — (5) — (6).  It  is  a  well  organised 
algorithm,  producing  the  same  result  as  previously  obtained. 

The  example  cited  above  has  the  property  of  being  discrete  which  is 
indispensable  in  the  dynamic  programming  algorithm.  Although  the  problems 
encountered  in  engineering  are  often  of  a  continuous  nature,  the  discrete 
interval  can  often  be  made  small  enough  to  approach  the  continuous  nature 
of  the  problem.  This  can  be  accomplished  with  the  use  of  modern  high  speed 
computers. 

The  Exploitation  Problem 

The  problem  of  exploiting  a  natural  resource  which  is  available  for  a 

total  of  N  years  will  be  considered  next  (9).  The  state  x  is  defined  as 

n 

the  amount  of  the  resource  remaining  for  processing  with  n  years  remaining. 
Suppose  that  at  each  year  a  decision  8  concerning  the  method  of  processing 
for  the  next  year  is  made  and  thus,  altogether'N  decisions  are  made.  The 
condition  of  the  resource  remaining  for  further  exploitation  and  the  profit 
made  during  exploitation  are  dependent  upon  the  decision  made  at  the  present 
stage.  For  example,  an  oil  reservoir  may  be  initially  water  flooded  to  in- 
crease the  rate  of  oil  removal  and  the  profit  rate;  but  this  early  flooding 
may  disturb  the  reservoir  and  in  the  long  run  reduce  the  total  amount  of 
oil  that  can  be  removed. 


Representing  the  sequence  of  years  in  the  order  shown  in  Figure  2,  the 
interval  profit  P  is  the  profit  made  during  the  n-th  period  and  is  a  function 
of  the  state  of  the  resource  at  the  beginning  of  the  n-th  stage  and  the  deci- 
sion made  at  that  stage,  i.e.,  P  ■  P  (x  _;8  ),  As  previously  mentioned,  the 
condition  of  the  resource  remaining  for  further  exploitation  depends  on  the 
decision  made  and  also  on  the  previous  ststa.  The  relationship  between  them 
may  be  denoted  by  T,  the  transformation  operator.  T  may  be  a  complicated  math- 
ematical equation  or  just  common  sense;  it  may  be  expressed  by 

x  ■  T(x   :6  ).  The  objective  function  is  the  summation  of  the  interval 
n     n+1  n 

profits 

The  value  of  the  objective  function  depends  on  the  initial  state,  x.  ,  of  the 
resource  and  a  sequence  of  N  decisions,  8  ,  8M  .,...,9  .  The  problem  is  to 
make  a  series  of  decisions  which  will  maximize  the  objective  function.  Such 
decisions  are  called  the  optimal  policies  denoted  by  0»,  Ojj_i>.>.»  8..  Thus 

o*C«ta>*  o(Xin;  e;;  o^,...,  e*)5 fH(^)  .  jtalj,  p(v  en)j   (3) 

f  (x  . )  is  the  total  profit  obtained  over  the  time  period  n  to  0  if  an 
optimal  policy  is  used  starting  at  stage  n.  Clearly,  if  we  find  ourselves 
at  a  state  x    at  a  certain  time  or  a  stage  n,  the  best  we  can  possibly  do 
is  to  follow  the  optimal  policy  from  the  state  x    to  the  end  and  thus  ob- 
tain a  profit,  f  (x  ,},  over  that  time  period.  Applying  the  principle  of 
optimality  the  dynamic  programming  algorithm  may  be  written  as 

£    (x  .  )  -  Max  {  P  (x  x.  ;  8  )  +  f  ,  (x  )  I  (M 

n  n+1    e   l  n  n+1'  n    n-1  n  |  w' 
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x  =  T(x   ;  6  )  (<) 

n     n+1 '     r,  w 


n  n+1'  n 


If  M  different  decisions  are  possible  each  year,  the  optimal  policy  and  ob- 
jective function,  equation  (3),  may  be  reached  by  the  repeated  and  simultane- 
ous use  of  equations  (1*)  and  (5).  The  procedure  starting  from  the  last  stage 
is  illustrated  as  follows: 


8i 


First,  for  a  selected  value  of  x  ,   the  values  of  P  (x ;  0  )  are  computed  for 
M  possible  8^  which  maximizes  P-^Xgj  0^)  is  the  optimal  decision  8*  corre- 
sponding to  the  particular  x    selected.     The  corresponding  x     can  be  cal- 
culated by  equation  (7).     Then  the  same  calculation  is  repeated  for  all 
possible  x2.     Eventually  the  so-called  dynamic  programming  table  may  be  con- 
structed as  shown  in  Table  1A.     Proceeding  to  the  second  stage,  equations   (U) 
and  (5)  are  written  as 

f2(x3)-Max|p2(x3;  e2)+fi(x2)  /  (6) 

x2  -  T(x3;  82)  (?) 

Again  assigning  a  value  to  x^   the  values  of  ^(Xjj  8  )  and  x  are  computed 
from  equations  (8)  and  (9)  respectively.  The  value  of  f  (x  )  in  Table  1A 
corresponding  to  the  computed  xg  is  added  to  P^x  ;  0  ).  The  optimum  profit, 
fgCxj),  for  the  particular  x  is  the  maximum  among  all  the  computed  values 
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Table      I 

Illustration     of     dynamic      programming 

algorithm 


Table    I. A  Table    LB 


x2    9*   f,cx2)  x,      x3   e*   f2(x3)  x2 


Table    l.n-1  Table    l.n 

xn     en-i  ^n-i<Xn>  Xn.,     Xnf,.en   fn(Xn+|)    X„ 


13 

of  JP  (x,;  6  )  +  £.(x?)l  for  all  possible  9  ,     The  succeeding  9?  and  the  x 
correspondent  to  the  computed  x„  are  located  in  Table  1A.  Table  IB  is  like- 
wise constructed  by  repeating  the  above  calculation  for  each  possible  x,. 
In  a  similar  fashion  all  N  tables  may  be  constructed  sequentially.  In  general, 

an  n-th  stage  table  has  the  form  of  Table  In.  With  a  particular  x..  .  "X. 
L  N+l    in 

given,  we  proceed  to  Table  N  and  find  0  and  x^  from  it.  Proceeding  to  the 
(N-l)  stcge,  Om  J  and  xl.   corresponding  to  the  initial  state,  x  are  ob- 
tained, and  this  process  can  be  continued  until  the  entire  sequence  of  the 
optimal  decision,  Gk»  8jj'_]_, ..  .,6£,  is  extracted.  This  algorithm  not  only 

provides  the  answer  to  the  particular  x    ■  x  ,  but  it  also  yields  a 

N+l    i^n 

series  of  solutions  for  a  whole  range  of  x  ,  meaning  that  the  dynamic 
programming  tables  contain  solutions  for  M  x  N  number  of  optimization  prob- 
lems. Knowing  this  is  especially  useful  in  chemical  engineering  problems. 
While  the  feed  condition  can  hardly  be  kept  at  a  steady  state,  we  can  de- 
termine the  effect  of  putting  on  or  taking  off  a  few  stages.  Furthermore, 
the  dynamic  programming  formulation  not  only  allows  a  problem  with  high 
dimensionality  to  be  replaced  by  a  sequence  of  problems  with  lower  dimen- 
sionality, but  it  also  provides  a  straight  forward  computing  algorithm.  The 
device  by  which  computational  solution  is  brought  about  is  through  a  system 
of  recursive  equations,  which  is  well  suited  to  high  speed  computer  machine 
solution. 

Another  way  of  ordering  the  series  for  the  stages  for  the  problem  just 
considered  is  presented  in  Figure  3.  The  algorithm  is  essentially  the  same 
as  before,  except  that  the  direction  is  reversed  in  numbering  the  states 
and  control  variables.  According  to  this  scheme, 
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OS  fN(xQ)  =  Max      PN(xo;  ex)  +  f^)  (Ua) 

While  the  former  way  of  numbering  the  system  is  more  canvenient  for  the  purpose 
of  a  computer  programming,  each  has  its  advantage.     The  latter  enables  us  to 
compare  dynamic  programming  with  the  maximum  principle. 

The  dynamic  programming  algorithm,  as  has  been  shown,  is  convenient  and 
powerful,  but  it  is  not  a  rote  to  be  thoughtlessly  applied  to  every  problem. 
One  of  its  best  features  is  that  it  permits  exercise  of  ingenuity  in  simpli- 
fying the  complexities  of  a  particular  problem  before  applying  the  basic 
algorithm.     Two  schemes  for  reducing  the  complexities  will  be  introduced  by 
means  of  concrete  examples. 

Lagrange  Multiplier 

Quite  often  the  use  of  the  Lagrange  Multiplier  can  reduce  the  dimen- 
sionality of  a  problem  and  enable  the  problem  to  be  handled  within  the  capa- 
city of  computers   (17).     Usually  the  Lagrange  Multipliers  are  applied  when 
there  are  inequality  constraints  or  restrictions  on  the  variables   (17). 
There  are  two  important  rules  in  applying  this  technique   (17):     1.  The  number 
of  Lagrange  Multipliers  to  be  introduced  must  be  equal  to  the  number  of  con- 
straining equations  and  2.  The  Lagrange  expression  must  be  equal  to  the  ob- 
jective function  plus  the  product  of  the  Lagrange  Multipliers  and  ',...'        is 
the  constraints.     These  constraints  must  be  in  the  form  of  an  equation  which 
is  equal  to  zero. 

Suppose  that  the  function  f  (x,y,z)  »  x  +  2y2  +  z2,  subject  to  the  con- 
straint that  x  +  y  +  z  -  1  =  0,   is   to  be  minimized.      The  Lagrange  expression 

(17)  is: 

w(x,y,z)  =  x  +  2y    +  z     +A(x  +  y+z-l) 

Where^A  is  the  Lagrange  Multiplier.  ConsiderX  as  a  new  variable. 
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aw 


1  + 


^-o 


#2  =  hy  +  A  -  o 
ay 


2z  +  A 


aw 

£?  =  x+y  +  z-l  =  0 
Solving  these  four  equations  simultaneously, 

A  -  -  l,  x  -  0.25,  y  =  0.25,  z  =  0.5 

2         2  2 

Thus,   the  minimum  value  of  f(x,y,s)B  x  +  2y    +2     -  (0.25)  +  2(0.25)     + 

2 

(0.5)  ■  0.625,  under  the  constraint,  x  +  y  +  z-l  =  0. 

Now  the  case  with  two  types  of  resources  represented  by  quantities 
x  and  y  respectively  will  be  considered  (k).     The  quantities  of  these  re- 
sources allocated  to  the  n-th  activity  are  designated  by  x  and  y 
respectively.  If  there  exists  an  interval  profit  P  (x  ,  y  )  which  is  the 

return  from  the  n-th  activity  due  to  respective  allocations  of  x  and  y  , 

n     n 

then  the  optimal  problem  is  to  maximize  the  objective  function  of  2N  vari- 
ables, i.e.,  the  proper  allocation  of  the  resources  for  an  efficient  utili- 
zation. In  other  words,  it  is  to  maximize 

0(xx,  x2,...,  xN;  yx,  %,...,  rj   -  g  \(\,   7n>  (10) 


subject  to  the  constraints 

"S7  x     ■  X 
J-*     n 
n-1 

x  ^  0 
n 

N 

si  y    •  y 

jrv  o 

n 

(11) 


17 


This  is  a  two-dimensional  allocation  process,  independent  of  the  value  of  N. 
Applying  the  dynamic  programming  algorithm, 


tfW)«  fN(x,y)  -  Max    P^,  X^...,^,  y^  y,,,..., 


V 


x     and  y     are  values  to  be  determined, 
n  "n 

For  N  -  1 

f1(x,y)   ■  P1(x,y) 

For  K22,   the  following  recurrence  relation  is  obtained, 

f  (x,y)  -  Qjto  ojx^  [Pg^jr,)  -  y»  -  V  y  -  yN)J 

The  solution  of  this  recurrence  equation  is  quite  involved  and  demands  a  large 
computer  capacity.  In  the  following  treatment  it  will  be  shown  how  the 
Lagrange  Multiplier  can  simplify  the  computation. 

The  following  equation  with  constraints  on  x  and  y 

o(xi'--v  yi v  ■  vvv +  vvv +  - +  y  w 

Xl  +  X2+  •"  +VX     •     \2° 
is  replaced  by  the  problem  of  maximizing  the  modified  function 

p1c«lf  tj)  *  p2u2,  y2)  *  ...y*j,,  yK)  -  ^i  *  y2  +  ••• +  yN_] 

subject  to  the  constraints 

x1  +  x2  ♦  ...  ♦  xjj  -  x     ,     xn20  (13) 


IS 


For  the  moment  consider  /\  as  a  fixed  parameter  and  assume  that 

p  U  i  y  )/y  -=*  o  as  y  — *■  °° 

n  n   n   n         n 
The  maximization  over  y  can  be  done  independently  of  the  maximization  over 

x  (It).  Thus 
n 

hn(V^)  =  hn(xn}  =  ^[VV  V  "^nl 
Tne  problem  is  transformed  to  that  .of  maximizing  the  function 

W  +  h2(x2)  +   •"+hN(XN) 

subject  to  the  constraints  of  equations  (13)  and  (111).     Repeating  the  pro- 
cedure by  varying  X  until  it  satisfies  the  restriction 

2  yn(A)  -  r 

n»l 

it  may  be  found  that  the  corresponding  optimal  policy  is  the  solution  for 
the  problem. 

Geometrically,  the  Lagrange  Multiplier  can  be  interpreted  as  the  slope 
of  the  searching  lines  in  the  complex  plane  (1*).  As  will  be  shown  in  the 
first  example  in  the  application  section,  the  Lagrange  Multiplier  can  also 
denote  price  or  relative  cost. 

Batio  of  State  Variables 

The  term  "homogeneous  nature"  implies  that  the  transformation 
functions  of  the  problem  are  linear  with  respect  to  each  other  (6).  This 
linear  property  simplifies  the  optimization  problem  in  that  the  ratio  of 
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two  state  variables  may  be  used  as  a  new  state  .variable,  and  thus   the  number 
of  state  variables  is  reduced  (6). 

As  an  example,   the  system  of  a  first  order  consecutive  reaction 
A-^B-^C,   will  be  considered. 

Hef erring  to  Figure  2  the  concentrations  of  A  and  B  which  are  denoted 

by  x  and  y  respectively,  are  the  state  variable  and  the  temperature  is  the 

controlling  variable.     The  n-th  reactor  has  the  volume  V     and  the  volu- 

n 

metric  flow  rate  of  the  stream  v  is  constant  throughout  the  multistage 
system.     It  will  be  assumed  that  there  is  no  density  change  in  the  process, 
which  is  essentially  true  for  the  liquid  phase  reactions.     The  holding  time 

is  e   .__» 

n         V 

Since  the  reaction  under  consideration  is  a  first  order  consecutive  reaction 
represented  by 

k(t)       h(t) 
A" >B >Q 


and  the  kinetic  equations  for  rate  of  appearance  of  A  and  B 


rA  -  -  k(t)x 

p  -  k(t)x  -  h(t)y 

the  material  balance  around  the  n-th  stage  for  each  species  is: 

vx    -  vx   +  (-r  )  V 
n+1     n    v  A   n 


x    ■  x  +  k(t  )  x  0  del 

n+1    n     n   n  n  ^AS' 
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and 


Vn*l  =  ^n  +   (-rB}  \ 


yn+l=yn     +  h^V  yn  "  k(V  Xn  8n  ^ 

The  production  of  y  for  a  definite  feed  condition  is  to  be  maximized  by- 
choosing  the  optimal  series  of  decisions,  6    and  t  .     It  can  be  seen  that 

n  n 

in  order  to  describe  the  system  completely  the  composition,  x  and  y,  must 
be  individually  specified  at  each  stage.     This  is  a  two-dimensional  prob- 
lem which  not  only  requires  complex  computational  scheme  but  also  demands 
a  great  deal  of  computer  memories.     However,  as  we  examine  closely,  we  can 
see  that  because  of  the  linear  transformation  equations  involved,   the 
problem  can  be  simplified  by  definirrg  a  new  variable 
x 


(17) 


Dividing  equation   (16)  by  equation  (15) 
**♦!  Xn  +  Xn  8n  k(V 


°n+l 


e     +  8     f  e    h(t  )  -  k(t  i) 
n        n  j    n        n'         ^  VJ 


1+0    k(t  ) 
n        n 


e         ■  e 


i  ♦  en  h(tn) 


n+l        n  (  1  +  6    k(t  ) 


0    k(t  ) 
n        a' 

1  +  8     k(t   ) 
n        n 


e         +  SO>)     JIV^V 

11+1       1  +  gnk(tnM       1*0    h(t) 


n+l  f1  +  6n  k(  VJ  +  8n  k(V  j  /  fr  *  °n  ^J 


(18) 
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This   equation  gives  a  functional  relation  for  determining  the  output  ratio 
at  any  stage  in  terms  of  the  feed  ratio  e^.     From  equations   (15)  and   (16) 


3 

n 


Vl  \  +  \  k<V  °n 

i  ♦  en  k(tn) 

Since  the  formation  of  B  is  to  be  maximized,   the'  objective  function  is: 


yi  "  yn+l  *   \    frn  "  yn+1)  (20) 

and  the  optimum  function,  f   (x       ;  y       ).   is 

M 

Because  of  the  linear  nature  of   the  equations,   this  must  be  proportional  to 
XjJ+1  and  depends  only  on  the  ratio  jr^j/J^  •  e^  i.e. 

VVl'  7w)*  XN+1  gN(eN+l}  (22) 

Now  applying  the  dynamic  programing  algorithm  to  the  function  f 

N 

N'  N  V  ) 

Combing  equations  (15),  (19),  (22)  and  (23) 

g  (e   )  »  Max   yN  "  yN+l    1  ,     , } 
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g   (3       )  -    Max  fe  k(tn}  -  8n  h\K  +  1 ,     J    ,  ,  , 

N  Vl   Cl\ — 1  +  enk(v~  ^v^y  ^^j (2U) 

This  is  an  equation  in  e  only,  and  can  be  solved  recursively  in  conjunction 
with  equation  (18). 

DISCRETE  MAXIMUM  PRINCIPLE 

In  this  section,  the  general  maximum  principle  algorithm  for  discrete 
systems  with  recycle  is  presented  first  ftj).     It  is  then  reduced  to  a  special 
form  with  no  recycle. 

A  series  of  multi-stage  ideal  backmix  reactors  with  recycle  from  the 
last  stage  to  the  first  stage  will  be  considered.     They  are  schematically 

shown  in  Figure  h.     If  the  superscript,  n,   indicates  the  stage  number, 

n.        cio 

xi  s,  i  »  l,   2,...,s,  are  the  state  variables  leaving  the  n-th  stage  and 

er>  r  =  1j   2>-">t,  are  the  control  variables  at  the  n-th  stage.     If  T.   is 
the  transformation  function  for  each  different  i,   then 

1-1,   2,...,s 

n  =  1,   2,...,N 
where  T?(x.       .   Onl  Is  a  >h««t  h=nri  f„~™  <•-»  *  /•_n"1.   .n-1  __n-l 


ira  Ii(xk"  I  V  is  a  short  hand  f°™for  I.(XJ     •  x?"1, . . . ,xn_1; 

The  feed  enters  the  system  at  a  rate  q,  while  the  recycle  stream  is  fed 
back  at  a  rate  r.     The  mixing  conditions  are  represented  by 

Xi  =  K{xl->  \)  >  i-l,   2,...,s  (26) 

where  x£  are  the  state  variables  in  the  feed. 


I. 
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The  problem  is  to  find  a  sequence  of  control  decisions  \en  {  ,  with  or 
without  the  constraints,   c("^9"— /";  n  =  1,  2,...,Nj  r  -  1,  2,...,t, 

which  will  maximize  x  . 
m 

The  essence  of  the  discrete  algorithm  for  the  maximum  principle   (5)  is 

to  introduce  a  new  set  of  variables  defined  by 
n,    n-1     n. 

i  -  1,   2,...,S 

n  =  1,   2,...,N 

The  Hamiltonian  H     is  formed  as  follows: 


(29) 


It  can  be  shown  that  (5)   (15) 
azi 

Then  the  optimization  problem  becomes  that  of  finding  a  seo.uence  of  ■[  Qn[ 

to  satisfy  the  following  conditions   (.$,15) 

s  n 

a     "  ZT      Sj   TjCxjj"   ;   e")  -  Maximum  (31) 

with  the  boundary  conditions   (15) 

xi  =  VV  V  (32) 
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Zi  "  K  Z3  J     \  "    ^  (33) 

1     Ji   J         a  4 

n  ■  1,  2,...,N 

i  ■  1,  2,...,x 


(1      i  -  ra 
S.     - 


where   &.     is   the  Kroneclcer  delta,  i.e.      ST      -i 


"1      i  =  m 
-0      i  /  m 

If  the  minimizing  sequence  instead  of  the  maximizing  sequence  of  control 
actions  is  to  be  decided,  the  basic  algorithm  remains  unchanged  except  that 
the  equation  H11  ■  max  is  replaced  by  Hn  »  min. 

Quite  often  the  maxima  sought  in  equation   (31)  are  found  to  lie  at 
stationary  positions.     This  is  equivalent  to 

n.  n-1      n. 

er  ;  2>    ■        J  ^  p     r     .  o  (3U) 

r     d-i     3  aoR 

r 
If  the  9's  are  under  the  constraints 

then,  for  each  n,   the  maximum  required  by  equation  (31)  may  be  found  at  some 
stationary  place  between  s/n  and   a  "  or  at  one  of  the  end  points,  (jn     or 
A    .     In  this  situation,   the  following  equations  must  be  used: 


r 


«r      ■     4 
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or   :   0     = 


J> 


This  algorithm  can  seldom  be  solved  analytically  (15).     Even  in  simple 
cases  considerable  ingenuity  is  needed.     Therefore,   only  a  very  general  work- 
ing procedure  for  finding  a  numerical  answer  to  a  specific  problem  according 
to  the  algorithm  is  presented: 

1.  Propose  a  sequence  6     ,   9  , ... 

2.  Propose  a  set  of  x'. .         i  =  1,  2,...,s 

f 

3.  With  a  given  x.,  solve  equations   (26)  and  (25)  forward  and  get  a 

new  set  of  x  .     Kepeat  the  calculating  using  this  new  set  of  x 

N  * 

until  the  subsequent  x.   are  close  enough. 

h.     With  these  x    and  011,  solve  the  equations   (27)  and  (33)  in  the  z's 
i  r 

backward  from  n  ■  N  to  n  ■  1. 

5.  With  these  x's  and  z's  compute  a  sequence  of  8n  from  equation   (31) 
and  return  to  step  2.  < 

6.  The  iteration  procedure  is   terminated  when  the  successive  computed 

N 
values  of  x    are  sufficiently  close,  possibly  together  with  the  con- 
m 

dition  that  successive  computed  strings  of  0  are  also  close  enough. 
Quite  often  a  modification  in  practice  may  save  considerable  work. 
Instead  of  taking  0  computed  in  step  5  directly' back  to  step  2,  it  may 
sometimes  be  desirable  to  regard  the  new  0  as  a  mere  indication  of  direction 
in  8-space  and  apply  the  methods  appropriate  to  steepest  descent  calcula- 
tions to  choose  the  best  length  of  step  in  this  direction,  and  possibly  even 
to  modify  the  direction  in  the  light  of  the  iterative  history. 
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For  the  process  without  product  recycle,  i.e.,  r  -  0,  equation  (32) 
reduces  to 


0   f 
x  -  x 

i    i 


1-1,  2,...,s 

and  the  basic  algorithm  reduces  to 


n 
[     - 

i 

2j£ 

_   n 
Jz. 

X 

n 

an" 

i 

af4 

n  =  1,  2,...,N     r  -  1,  2,...,t     1-1,  2,...,s 

0^  are  determined  where  Hn  -  max,  with  the  boundary  conditions 


n 


0         f 
x.    -  X. 

i         i 

i     \0, 

i"m 

i^m 

(36) 


(27) 


n-l       A,    9*K*i  gj       n 
z.       =   21  " — =-* — — —  z^ 

1         3-1        ^x"-1  J 

n        A        n     n,  n-l       n. 


(29) 


(30) 


(36) 
(37) 

Hie  general  working  procedure  is  the  same  as  stated  for  the  case  with  product 
recycle  except  that  the  step  of  guessing  a  set  of  x1?  may  be  omitted  saving  con- 
siderable labor  in  matching  the  JT  by  trial  and  error. 
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It  is  worthwhile  to  note  that  for  each  assigned  value  of  x  ,  the  corre- 
sponding calculated  values  of  x",  n  =  0,   1,  2,...,N-1  are  the  optimal  state 
variables  corresponding  to  the  initial  condition  x    computed  in  each  run  of 
trial  calculation. 

If  the  final  values  of  certain  x's  other  than  x    are  specified,   the 
boundary  conditions  of  tquation  (37)  are  changed  tot 

N       [  1,     i=m 
V[0,    Wrf.p  W 

xp  "  b  (39) 

Otherwise  the  algorithm  is  the  same  as  before. 

APPLICATIONS 

Cross-current  Elutriation  of  Sludge 

The  purpose  of  elutriation  in  the  sludge  treatment  is  to  reduce  the 

alkalinity  of  the  sludge  with  wash  water  or  wash  solution  (18).     Here,   the 

term  "sludge"  indicates  the  sludge  in  wet  slurry  form  after  the  filtration 

process.     The  alkalinity  of  the  sludge  is  usually  expressed  as  mg/liter  of 

calcium  carbonate  equivalent  to  the  determined  alkalinity.     The  elutriation 

can  be  carried  out  practically  by  repeated  cross-current  washings  as  shown 

in  Figure  5  (18).     The  sludge  is  fed  at  a  constant  rate  corresponding  to  liter 

of  water  per  unit  time  and  with  alkalinity  E     mg/liter.     In  the  tank,   it  is 

o 

kept  in  suspension  by  air  or  mechanical  agitation.     After  thorough  washing, 
the  sludge  and  wash  water  are.  separated.     The  out  wet  sludge  solids  and  wash 
water  should  have  the  same  alkalinity.     The  outlet  sludge  solids  become  the 
feed  to  the  next  tank.     The  wash  water  is  available  at  the  alkalinity  w  mg/liter. 
It  is  fed  into  the  tank  at  the  rate  of  A  liter  per  unit  time. 
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For  a  system  with  a  finite  number  of  stages  and  a  definite  input  alkalinity 
E  ,  two  types  of  problems  may  arise: 

1.  With  a  fixed  total  amount  of  wash  water  & ,   how  can  this  amount  of  wash 
water  be  distributed  amoung  the  stages  in  order  to  maximize  the  reduction  in 
alkalinity,  E  -E  ,  of  the  sludge? 

2.  If  the  cost  of  providing  and  pumping  the  wash  water  is  a  fraction, 
A   ,  of  the  value  to  be  gained  by  reducing  the  alkalinity  of  the  sludge,  how 

should  the  flow  rate  of  wash  water  to  each  stage  be  regulated  to  maximize  the 

net  profit,  E  -  E  -Aq  ? 
'  0   N    N 

It  will  be  shown  later  that  the  solution  of  the  second  problem  also 
yields  the  solution  of  the  first  and  is,  in  fact,  the  easiest  method  of 
solving  the  problem.  The  realtive  cost,  A,  plays  the  role  of  the  Lagrange 
multiplier  (U,7). 

It  will  be  shown  below  that  dynamic  programming  and  the  maximum  principle 
algorithms  yield  the  same  solution  for  this  particular  problem.  The  optimal 
decision,  as  will  be  shown  later,  is  to  distribute  the  wash  water  evenly  among 
all  the  stages. 

For  the  convenience  of  formulating  the  problem,  the  following  notations 
are  defined: 


(1»0) 
E=  x 

It  can  be  recognized  that 

e  =  * 
q 

The  physical  condition  demands   that 

a^>  x,  j> . .  .> x  ^. .  .'>  x  ,>  w 

—    1—        —    n—  H~ 
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Dynamic  Programming  Solution.     The  material  balance  around  the  n-th  stage 
with  respect  to  alkalinity  yields 

qE     ,   +  a  W     =     RE     +qE 
n-1        n  n  n      ^     n 

n-1        n  n  n        n 

x    ,       9    w 

n 

This  represents  the  transformation  function,   T(x  :   S  ). 

'        n'     a 

Interval  profit  per  unit  volume  of  the  wet  sludge  at  the  n-th  stage  is, 

P     =  x        -x     -/ie 

n        n-1        n  n 

Therefore,  according  to  equation   (l;a)  and  Figure  3, 

*  (  \      N 

0   B  f   (a)  =  Max    )  a  -  x     -  X  "ST     9 

N  f9nl   '  ^      n 


Max 


\ )!  ^--.--v] 


Max      \  P  (a,  0  )  +  f       (x  ) 
<    »v   '     r         N-lv  l' 


f   (a)  '  Max     J    a  -  x     -  \q    +  f       (x  )  , 
N  ^  11         n-1     1    I 


(1*2) 


Clearly,  if  there  is  no  stage  existing,   the  profit  will  be  zero,  i.e. 

f0M  -  0  (1,3) 

For  the  one-stage  system,  shown  in  Fig.  6A,  equation   (1|2)  is  reduced  to, 


Va)"*f  y-\-\  +fo(vf 
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Combining  this  with  equations   (li3)  and  (Ul) 

a  +  0_  w 

Olli) 


fl(a)-Max       a-       _r_l_     .  ^ 


By  straightforward  differentiation,  it  yields  the  optimal  decision  0  as 


o   =  Ctsi)     -  i  (W) 

1         /\ 


Inserting  equation   (US)  into  equations    (Ul)  and   (UU) 


x     ■ 

1 


A(a  -  w) J     +  w  t  (U6) 

1  1 

f1(a)  -  a  -  w  +  j\-  2  ^2(a  -  w)2  (U7) 


For  the  optimal  two-stage  policy,  referring  to  Figure  6B  and  employing 
equation  (U2), 

f,(a)  -  Max  J  a  -  x,     -  A  0     +  f ,  (x  )  I  (U8) 

Combining  this  equation  with  equations   (1;1)  and  (li7) 


f2U) 


Max    j    a-i     -  Ae       +x       +A-w-2/i(x     -  w)        | 


f2(a)  =  Max    S    a  + /i  -  w  -Ao, 

By  direct  differentiation  the  optimum  point  is  found  as 
7 

(5D 


xx  4A(a-w)2]j  +  w 


3b 


f2(a)  -  a  -  w  +  2  A  _  3  /\3     (a  -  w)  3  (S2) 

According  to  the  principle  of  optimality  the  second  tank  should  use  the 
optimal  one-stage  policy  with  respect  to  its  input  x  .  Therefore,  from 
equation   (it5) 

x-w    * 
\  =   (7T)         -  1  (Wa) 

Inserting  equation  (Si)  into  the  equation  above 
1 

V^-1  (S3) 

Comparing  this  with  equation  (SO) 

e2  =  9l 
Thus  for  the  two-stage  system,   it  has  been  proved  that  the  optimal  policy  is 
to  divide   the  washwater  equally  between  the   two  stages. 

In  general  for  the  N-stage  system  it  can  be  deduced  that 

fM(a)  =  a-w  +  NA-(N+l)AI:fI(a.w)^  (&) 

N      (~j-)  -1  (55) 

8=8  -,      n 

n        n-1  n  =  1,  2,   ...,  n 

The  optimal  policy  is  always  to  distribute  the  wash  water  evenly  among  all 
the  stages. 

Maximum  Principle  Solution.  This  problem  is  a  simple  case  without  recycle. 
Following  the  notation  of  the  maximum  principle  algorithm  and  from  eauation 
(Ul) 
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X     1   K  1  + 


"-1  +  en  w 


? 


Defining 


n   „  ,  n-1   n^    n-1    n-1    n   \  n 
x  -  T  x   :  8  )  =  x,   +  x    -  x"  -  A  8 
2    2  tc         1     1     1 


(56) 


(57) 


.N 


H 


n-1    n 


it  can  be  shown  that  x"  =  2  (x    ~  xi  ~  ■"*  J  ^  *^e  objective  function 

2   n=l   !     x         M 

to  be  maximized.  The  problem  is  then  to  maximize  x  with  the  boundary- 
conditions 


1 

0 
x  =  0 
2 


From  equation  (27) 


n-1 


1  +  8 


1  +  8 


n  =  1,  2,...,N 


n-1    n 
z    =  z 
2      2 

From  equation  (37) 

N   „ 
z  -  0 

1 

z  -  1 
2 

From  equations  (59),  (60)  and  (61) 

n-1   n 
z    -  z  >=  1 
2      2 

n 
•z 

„n-l  . 


n-1,  2,...,N 


B 


1  +  0"     1  +  8" 
From  equations  (28),   (56)  and  (62) 


(58) 

(59) 
(60) 


(61) 

(62) 

(63) 


3o 


H0 


X?"1  +  9°  W 


n-1    n 
x,   +  9  v 


i (- — H  ♦  4(4  +  4"1  -  - — -n-  -  A  en) 

l  +  o     "  *    •■-      i  +  er 


h"  -  (a°  - 1)  (.1 ; 


♦  x*'1  ♦  x^-1 


A 


^0"         (i+  Qaf  '   /l  (64) 

Since  A  imposes  a  restriction  on  011  automatically  the  nature  of  the  physical 
situation  assures  the  existence  of  the  extremum  points  at  which  911  make 
equation  (61;)  equal  to  zero  (IS).  Thus 


(w  -  Xj_  )  (Zl  _  i) 


J 


A 


Solving  equation  (65)  for  z" 


(1  *  6n)  A 
n-1 


+  1 


(65) 


(66) 


For  the  (n+l)-th  stage 


1  n  -1- 


W  -  X. 


(66a) 


From  equation  (63) 


n+1 


1  i  +  en+1 


n+1 


l  +  e" 


Substituting  equations  (66)  and  (66a)  into  this  equation. 


(i  *  e") 

n=I 

w  -  x 

1 


(i  j  on+1) 


(67) 
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From  equation   (56) 


8    +  1  "  -h (68) 

x,   -  w 


n+1  __i 

8       +  1  °    x^+i  -  «  (68a) 

Combining  equations   (67),    (68)  and   (68a) 

(x^1  -  w)   (x^+1  -  w)  =  (*£  -  w)2  (69) 

The  general  solution  of  this  finite  difference  equation  is   (19): 

x^    =  c(d)n  +  w  ,       n  -  0,  1,...,N+1  (70) 

From  the  boundary  conditions,   the  constants,   c  and  d,   can  be  evaluated. 
From  equation  (58) 

c  ■  a  -  w 
Equation  (69)  then,  becomes 


(71) 


3C°  =  (a  -  v)  (d)n  +  w  (72) 

At  the  N-th  stage,  by  equations  (68)  and  (68a) 


H-1 


6M  +  1  = 

X    -  w 

N 
x  -  w 

1 

Combining 

equations 

i  (61)  and 

(66) 

N 

e  +  l  - 

H-1 
x,   -  w 

<— ; ) 

A 

Thus: 


(68b) 


(73) 
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(xj  -  w)2  =  (x*"1  -  w)  A  (7U) 

Inserting  equation  (72)  into  equation  (7li)  yields 

d  -  (If1  (75) 

a-w 

The  general  solution  for  x  is  then  obtained  by  substituting  equation  (75) 


into  equation  (72) 

n 
N+l 


2?  .  (a  _  w)  (-A_) 
]_  a  -  w' 


(76) 


Ins  er  tin;  " 


g  x.^  from  equation  (76)  into  equation  (56) 
1         1 


N+l 
0n  =  (S-^-2)    -  1  (77) 

N+T       ^^ 

xN  =  (A)     (a  _w)      +_w  (78) 

By  a  simple  summation  of  equation  (57) 

»J  ■«£♦»$ -aJ-sA*  (79) 

Combining  equation  (79)  with  equations   (58),    (76),    (78)  and  (77) 

H  1 

N  v  N+l  N+l 

x     -  a  -  w  +  N/\-  (N+l)   (/\)  (a-w)  (80) 

Since  8  is  independent  of  stage  number,  as  can  be  seem  from  equation  (77), 
the  optimal  policy  is  to  distribute  the  wash  water  evenly  among  all  the 
stages. 

For  this  simple  problem,  both  the  dynamic  programming  and  maximum 
principle  approaches  result  in  identical  analytic  solutions  although  the 
latter  approach  is  more  direct  than  the  former. 


3* 


These  solutions  are  for  the  second  type  of  .-the  problems  listed  at  the 

beginning  of   this  section,  but  it  also  contains  solutions  for  the  first   type, 

i.e.,   the  case  with  fixed  C3  .     Since  the  optimal  policies  for  different  values 

K 

of  A  ,   which  is  a  paramter,  give  Q  =  2±      °  =  Qn(A)>  bT  searching  through 

N       n=l  " 

different/A,  suchAcan  be  found  that  Q_(  /[)  matches  the  given  Q_  as  described 
in  the  previous  section  describing  the  use  of  the  Lagrange  multiplier. 

Step  Rocket  Problem 

A  step  rocket  is  a  multistage  rocket  system  employing  the  technique  of 
jettisoning  structural  weight  during  flight  to  increase  performance.     Figure 
7  is  a  sketch  of  a  three-stage  step  rocket  system  (20).     The  rocket  weight  is 
defined  as  the  initial  gross  weight  at  the  start  of  a  stage  which  is  the  period 
of  time  starting  just  after  the  jettison  of  one  structural  package  and  lasting 
until  just  after  the  jettison  of  the  next.     The  difference  of  weight  between 
two  consecutive  stages  is   called  the  step  weight,  and  weight  remaining  after 
the  jettison  of  the  last  stage,  the  rocket  payload. 

The  purpose  of  the  rocket  system  is  to  attain  a  specified  final  velo- 
city carrying  a  desired  rocket  payload   (21).     If  the  number  of  stages  and 
the  material  and  propellant  for  constructing  the  rocket  system  have  been 
decided,   the  problem  is  to  determine  the  optimum  size  of  each  stage  so  that 
the  gross  weight  of  the  system  is  minimized.     A  two-stage  optimization  problem 
has  been  solved  by  Goldsmith  (22)  by  using  calculus  and  a  three-stage  problem 
by  Schumann   (23).     To  make  hand  calculation  possible,   certain  simplifying 
assumptions  have  been  made  in  their  solutions.     For  many  problems  the  assump- 
tions are  quite  reasonable  and  the  results  are  found  to  compare  favorably, 
with  actual  practice   (20).     The  simplifying  assumptions,  however,  must  be 
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verified  by  comparison  with  solutions  obtained  by  removing  all  or  some  of  the 
assumptions.  Also,  many  practical  problems  may  have  to  be  solved  by  use  of 
other  approaches. 

Dyke   (20)  employed  the  dynamic  programming  algorithm  to  solve  the  prob- 
lem, which  will  be  summarized  below.     A  working  scheme  suitable  for  the 
machine  computation  of  this  algorithm  and  the  maximum  principle  solution  are 
also  presented  in  this  section.     These  approaches  not  only  can  be  used  to 
check  the  simplified  solutions  mentioned  before  but  should  also  yield  reliable 
results  for  rocket  systems  with  more  than  three-stages. 

Before  introducing  the  general  transformation  equation,   the  following 
terms  will  be  defined   (20): 

c  ■  propellant  exhaust  velocity  of  stage  n 

f  (V)  =  minimum  weight  of  rocket  n  achieving  velocity  V 

M  =  ratio  of  initial  thrust  to  weight  of  rocket  n 

n  ■  number  of  stages 

V  ■       velocity  ideally  added  by  the  rocket  n;  design  velocity  of 

rocket  n 

W  «       rocket  system  payload  weight 

Jj 

w  =  initial  gross  weight  of  step  n 

K^j  =  that  portion  of  stage  n  jettison  weight  which  is  constant 

(/nwn  =  portion  of  stage  n  jettison  weight  dependent  of  step  weight 

'-^nVSi  ■  portion  of  stage  n  jettison  weight  dependent  on  thrust 

G~n  n  »  total  jettison  weight  of  stage  n 

xn  =  initial  gross  weight  of  rocket  n 

8  «  velocity  ideally  added  during  stage  n 
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a  =       -Ah      -J 
a  ran       <*  n 

b  =       1  -</ 

n  '  n 


n 


d  "         W 

n  xn 


It  is  assumed  that  M_  is  constant  for  the  stage  n.     The  following  rela- 
tion holds  for  each  stage   (20). 

(fi      =    A     *      +   /5n  M.  x    +  w  (81) 

0  n  n  1  n    n        ra    »     n        xn 

i.e.  for  stage  n,   total  jettison  weight  is  the  sum  of  the  jettison 
weight  which  depends  on  the  step  weight  and  step  thrust  plus  a  constant  tern. 


This  assumption  is  quite  reasonable   (20).      d   ,    iff    and  w      are  con- 

(82) 


straints  for  n-th  stage.     From  the  definition  of  x 

n 


The  well-known  rocket  equation  yields   (20) 

Xn  Sn 

x     .   +  fa  wn       =  exP   <ol>  ■  (83) 

n-1     v  n    n  n 

Combining  equations   (8l),  and  (82)  and  (83) 
(1  -<}  Vl  +  Wxn 


%-£*»-<* 


exp(--4-  ^„Mn-^n 

n 
b    x     ,  +  d_ 
n    n-1        n 


(BU) 


exp  (  -  — )  +  a 
n  " 

Dynamic  Programming  Solution  (20).  Applying  equation  (h) 

f  (V  )  -  Min  J   n  *   an  n-1   n-r   (  (8£) 

n  n  oiSiv    )         ;     er 

a   [     exp  (  _  -ii)  +  a 


13 


v  ,  »  v   -  e 

n-l        n        n 


(86) 


For  the  last  stage 


w 


dl  *  ai  WL 

V, 
exp   ( i)  +  a 


(87) 


An  outline  of  the  computer  flow  chart  for  solving  these  recurrence  rela-  , 
tions  is  shown  in  Figure  8  and  a  Fortran  1620  program  is  given  in  Appendix 
A. 

Maximum  Principle  Solution.  Applying  the  maximum  principle  algorithm,  the 
problem  can  be  reformulated  as  shown  in  Figure  9  which  corresponds  to  Figure 


7.  From  equation  (81i) 

, n  n-l   ,n 
b  x,   +  d 


exp  (  -  0y)  ♦  an 


(8Ua) 


Defining 


2    2 


n-l   _n 


(88) 


the  problem  is  then  to  minimize  x  with  the  following  boundary  conditions 

x  »  W 

1    L 


0   „ 
x  »  V 

2 


*j[-0 


(89) 


I* 


Start 


Read  in  an,  b",c",dn,wLyN 


Set  V,  =  9| 

fltn   B 


Compute  f 


Set    A-»Vn(r.22) 


bo 

J~ 

J 

C 

i 

©n 

Vn-l  =  Vn-Qn 


No 


Compute  f 

i 

fn  <<M 

No 

,S      YPS 

fn  —  /. 

* 

en+A— en 

i 

Gr 


D  1 


Yes 


Punch  Vn,^8i 
correspondir.gj  6n.Vn-i 


Vn*A— =-Vn 


No 


Vn=-B  ? 


~v? 


n+  i  ■ 


-n 


n  =  n  ? 

A     Yfis 

Stop 

No 


Fig.  8.      Outline    of    computer    program    for 
rocket     problem. 
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It  can  be  shown  that 

N      n-1        n        N 
/^  (x„       -  x2)   =2  0 

n=l  n=l 


N 
n-l 


Applying  the  end  point  specification  case  of  the  maximum  principle  algorithm 
equations  (27),  (28),  (38)  and  (39)  become 


(90) 


Thus 


b     z 

n-l 

1 

z 

1 

exp(-  £.)  *  an 
cn 

n-l 
z2 

-*2 

*• 

i 

*2  = 

Z2       Z2 

n-l,  2,...,N 

/  n  n-l        n 
b  X,       +  d 

H11- 

znJ        X 

(91) 


(92) 


exp(  -  5_)  +  a«» 


Because  of  the  introduction  of  the  second  state  variable,  x  .   the 

'     2' 
constraints,  0i8nflV,  are  automatically  removed  as  indicated  by  the  rela- 

N 
tion     2T    e"  =  V  and  the  fact  O1^  0. 

The  optimal  sequence  of  (enf  therefore  must  occur  at  the  extreme  points. 
That  is: 


U7 


3  8*    2   =n 


n  (b^"1  -  dn)  exp(-  £) 


=  0 


Lexp(— ^)  +  an 


Solving  for  z 


V 


exP  I-  — )  +  a 
c 


1    ,  n  n-1    n,     .   9n 
(b  r.   +  dn)  exp  (  -  — ) 


(93) 


likewise 


n-1 


n-1 


n-1      n 

,  8 ,    n-1 

exp  (-  -jj^jJ  +  a 


n-1  n-2   ,n-l 


(b   x 


-a-1 


+  d   )  exp  (-  2_  ) 

1  cn"i; 


(A) 


Inserting  equations  (93)  and  (9U)  into  equation  (90) 


n-1 


0   i    n-1 


n-i     ,  w   . 
(exp(-  --j)  * 


n-1  n-2  .  ,n-l. 


rfl-1 


(b   x    +  d   )  exp( _) 

1  n-1 

c 

From  equation   (81;a) 


exp  (-  2L) 
c 


bV"1  -anx"  t  dn 


"5 

, n  n-l        n 
0n   ,   .     n        b  xl      +  d 


exp   (-  -— )  +  an 


bV     J  exp  (-  t)  ♦  a» 


(b^"1  +  dn)  exp   (-  £) 


(95) 


(96) 


(97) 


likewise 


U8 


n-l    n-2        n-1  n-1        n-1 
n-1        b\  -a      ^      +d 


ex*>  (-  stf 


n-1 
xl 


en-l  ,       b^-^d-1 

o 


n-1 


(98) 


(99) 


Combining  equation   {96),   (96),    (97),   (98)  and  (99) 


la"/  • 


For  the  N-th  stage,  from  equations  (90)  and  (91) 


(100) 


N-1 


,N 


exp(-  0  )  +  a 


Substituting  equation  (91)  into  equation  (93) 


exp(-  *,) 


V 


exp  (-  Sjj)  *  aK 


2  N  N-1         N 


(bJj      +  d   ) 


(101) 


(102) 


For  the  (N-l)-th  stage,  from  equation   (9h) 


N-1 


exo( ) 

N-l; 


N-1 


N-1 
/      0        .  N-1 

6XP(-  TZ}   +  a 


2  N-1     N-1  N-2         N-1 

z,      (b       x         +  d       ) 


(103) 


U9 


Dividing  equation   (102)  by  equation   (103)  and  inserting  equation  (101)  for 
N-1 


1 


/     oW,  iil,  N-l  n_2        N-l  /-I 

exp(-~)  bcN(b      ^      +  d       )exp(-°      ) 


c 


exp(-^)  ♦  aK 
c" 


(10U) 


c  -J  l 


The  working  scheme  for  digital  computation  is  as  follows: 

1.  Since  x^  and  x°  are  known,  x  and  x^  can  be  evaluated  from  equa- 
tions (61ia)  and  (88),  by  assuming  a  value  for  0  . 

2.  Next,  x"  can  be  calculated  directly  from  equation  (100)  for  stages 
2  to  (N-l).  Then  0n  and  x£  are  calculated  from  equation  (8Ua)  and  equa- 
tion (88)  respectively, 

3.  Finally,  Q  is  obtained  from  equation  (10U)  for  the  last  stage 

N 
(N-th  stage).  If  x'2   calculated  by  substituting  the  computed  8N  into 

equation  (68)  is  not  zero,  the  procedure  is  repeated  from  the  first  step 

by  assuming  a  new  value  for  9^. 

i  N 

k.     If  the  x'    thus  obtained  is  zero,   this  sequence  of  911  is  the 

optimal  control  actions  sought,     x'1  can  be  secured  from  equation  (8Ua) 

A  similar  working  scheme  may  be  followed  by  assuming  a  value  of  x1 

instead  of  9  . 

Stagewise  Biochemical  Reactors 

In  this  section  a  generalized  optimization  problem  of  first  order 
consecutive  biochemical  reaction  is  considered.  It  is  numerically  solved 
in  detail.  A  comparative  study  of  the  two  approaches  is  also  made. 


SO 


In  many  of  the  biochemical  reactions,   the  catalyst  is  an  enzyme  of  pro- 
tein origin   (21;).     Although  the  catalytic  behavior  of  enzymes  is  quite  com- 
plicated,   the   two  most  influential  factors   on  the  activity  of  anzymes  are 
temperature  and  pH  value   (210.     An  increase  of  temperature  or  pH  value 
will  enhance  the  activity  of  the  enzyme;  but,  beyond  a  certain  extent  the 
enzyme  would  be  deactivated  by  further  increase  of  terrperature  or  pH  value. 
This  feature  of  the  biochemical  reaction  is   different  from  that  of  nonbiochem- 
ical  reactions  whose  specific  reaction  rates  change  according  to  the  Arrhenius 
rate  expression.     The  fact  that  there  exists  an  optimum  condition  for  the 
enzyme  activity  is  well-known  (21;, 25).     Though  the  exact  functions  of  the 
activity  versus   temperature  and  pH  value  should  be  acquired  under   the  actual 
operating  condition,  the  general  form  can  be  represented  by  the  unimodal  graph 
similar  to  that  shown  in  Figure  10  which  is  a  sketch  of  the  following  hypo- 
thetical functions   (25): 

k(0)  =  -  3S2  +  10  9  +  12  (105) 

h(0)  =  -  Q2  ♦  98  -  6  (lo6) 

where  Q  is  the  variable  representing  either  temperature  or  pH  value,  and 
h(8)  and  k(9)  are  the  activities  of  enzymes.  ' 

The  following  consecutive  biochemical  reaction  system  will  be  con- 


sidered 


A 


k(9) 


■j>b  Mei_^c 


Enzyme-!  '  Enzyme  2 


This  reaction  system  is  fairly  common  in  the  fermentation  processes.     Some 
of  the  examples  are  listed  in  Table  2    (25).      The  enzyme  indicated  is   the 
catalyst  for  each  reaction.     Many  of  them  are  of  commerical  importance. 
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©■ 


Fig.  10.     Graphs     of 


5  6  7 

k(e)    and    h(e) 
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The  raw  material  A  is  expensive,  B  is  waste  and  C  is  the  de- 
sired product.   If  a  series  of  ideal  equal  volume  backmix  reactors 
shown  in  Figure  2  is  available,  the  problem  is  to  find  an  optimal 
decision  sequence  <  0  >   in  order  to  maximize  the  production  of 
C  for  a  fixed  input  A. 

The  feed  rate  is  maintained  constant  and  thus  the  residence 
time  in  each  stage  is  equal,  let  it  be  t .  A  solution  is  presented  in 
a  general  form  first  and  then  a  three-stage  optimization  problem  is 
numerically  solved  for  the  case  with  the  residence  time  in -each  stage 
t  =  0.01. 

Dynamic  Programming  Solution.   Let  x,  y  and  w  represent  the  con- 
centrations of  A,  B  and  C  respectively.  The  material  balance  for 
the  n-th  stage  is 

xn+l  "  xn  +  k(6n>  V  (10?) 

^n+1  "  ^n  +  t  [h^J   ^n  '  k(8n>  xn]  '(10S) 

w„^  ■  w.  -   »(°J  *„*  (109) 


Defining 


y 


n+l 

n+l   x  . 
n+l 


w   1 
n+1   xn+l 


(110) 
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TABLE  3 

Some  examples  of  consecutive  biochemical  reactions  (25; 

,    _  ,.  .        Cellulase    _  . .  ,  .  Cellobiase     „, 

1.  Celluloses  »*    Cellobiose  ■»  Glucose 

2.  ^lal*°?e       Trehalase     >   Glucose     Zvmase Alcohol  +  00, 

Trehalose                 -^  2 

-i        t    u          Inulase    _     ..   ,  Zymase    _   ,,   ,  ,    „. 

3.  Insulin     S»-    Fructose    ■ &■     Alcohol  +  CO, 

2 

Albumoses      _    .  „               Hydroxy-acid 

i        n  Erecsin         ,  .    ...   Desaminases  , J  J 

4.  Peptones         '  «»"    Ammo  Acid  — ^ £*>     + 

Peptides  Ammonia 

v   ...        Xanthoxides  _      ...   ,  .  .     Uricase ...   .  . 

5.  Xantnme    »-    Uric  Acid   - ^  Allantom 

T               Fructose  _ 

t        o  Invertase  Zymase   ,.   ,  .    _, 

6.  Sucrose     w        +       » »■  Alcohol  +  CO, 

Glucose 


5<+ 


and  dividing  equations  (108)  and  (109)  by  equation  (107)  and  com- 
bining the  resulting  equations  with  equation  (110) 


u   fl  +  h(6  )t]   -  k(Q  )  t 
n  L v  n'  J v  n'  i,,,\ 

Jn+i  ■ 1  J  k(e  )t <m) 


II 


t  h(en)  un 


Vn+1 1  ;  k(0  )t  <112) 

v  n' 

i'rom  equation  (109)  r  x 

y  h(e  )t  f y„+i  +  k^e  >*  <i  ;  III  ;t> 

a-  h(e  )t 


*.     ,      X  .  "^V  "  x  .         1  +  h(0  )t 
n+1       n+1  n+1   /  *  n 

h(9  )t   (  k(S  )t 


Xn+1  T-Thtojt  jUn+l   +   1   J  k(Sn)t  <n3) 

The   objective   function   of  f„(x„T   .  .   y..  .  )    is 

N  N+1   N+1 

N 
VXN+1'  W  "  *\2j  <wn  "  Wn+1>  (U^ 

It  can  be  seen  from  equation  (113)  that  equation  (11^)  depends  only 
on  the  ratio  of  1^  .  and  y„  .  which  is  il  .<   This  reduction  in 
t'ne  number  of  state  variables  from  two  to  one  is  possible  because 
of  the  homogenous  nature  of  the  first  order  consecutive  reaction. 
Therefore: 


fN(xN+l'  yN+l}  "  *H+1  %(UN+l)  <U5) 
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Now,  applying  the  dynamic  programming  algorithm  to  the  function  f. 


H 


VXH+1«    yN+l^    "  M*X  f<wN   "  WN+1^    +   fH-l(V    V 


Combining  equations    (113),    (H5)    and   (116) 


(116) 


th(eK)-t  k(S  )t 

i  +  h(oH)t  <un+i  +  i  ;  k(oN)t)-+  i  +  k(eN)t  Su-i^V 

(117) 


The  general  recursive  equation  j. or  the  n-th  stage  is, 


/  h^Sn^  k^0n^  1 

sJun+i)  "  Kax  ti  +  "(e  )t  <un+i  +  i  ;  k(e  )0  +  i  J  k(e  )t  sn-iK) 


(118) 


Equation  (118)  contains  only  one  state  variable,  u,  and  it  can  be 

solved  recursively  in  conjunction  with  the  following  equation: 

u  .  |l  +  k(0  )t]  +  k(0  )t 
n+1  L v  ny  J ^  n/ 

n  "      1  +  h(Q  )t 


(119) 


A  flow  chart  for  solving  the  recurrence  equations  is  shown 
in  figure  11,  and  a  Fortran  1620  program  is  given  in  Appendix  B, 
The  results  for  u  ,    ranging  from  0.03  to  2  are  presented  in  Figure 
12.  About  12  hours  were  required  to  obtain  such  results  using  the 
increment  of  0,  0.1,  and  that  of  u,  0.03.   The  range  of  6  considered 
was  from  1  to  5.  The  computer  used  for  this  calculation  is  IBM 
1620  with  60,000  memories. 
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Start 

A 

Read   in  ,  T 

i 

Set  fn-,  =  0 
Xn_,     (AtoB) 

x 

Set 
A— xn 

A 

C 

em — ^-TnA 

i 

'  ■ 

Set 

C =-9n 

.1 

Compute 

Xn-..    Pn 

A 

Evaluate   fn 

j. 

fn^fnA? 

H 

IYes 

f  n       "~  fn  A 

:.. 

ml 

en+  a— en 

j. 

No 

e„^-D? 

IYes 

Output,  fnA,  Xn,  and 
corresponding  Gn,  Xn-i 

A 

Xn+  A— Xn 

I      . 

No 

Xn>E    ? 

- 

1  Yes 

n+i— n 

.1 

No 

n  ?tN? 

,1  Yes 

Stop 

Fig.  1 1.      Flow     chart     for    dynamic     Programming 

solution      of      bjocnemrcal      reaction     problem. 


5? 


0.03   0.18     a33  0.48a43    0.73  0.03   1.08    1.23  1.38    1.53     1.68    1.83    1-98 

Ratio       of      concentration 

Fig.  12.  A  graphical    representation   of    the    results 
of    the     dynamic     programming    solution 
for    biochemical     reaction. 
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Maximum  Principle  Solution.  The  notations  lor  the  maximum  principle 
algorithm  are  adopted,  x  ,  x  and  x  being  used  to  represent  the  con- 
centration of  x,  y  and  w  respectively. 

.from  equations  (107),  (103)  and  (109) 


n-l    n  .    ,  ■,  An-,.  ( 
k.       -   x.  / 1  +  k(9  )t 


n-1    n 


nN  n-1 


1  +  h(en)t     . -L-i-1 — 

J     1  +  k(9n)t 


x  "  =  x  -  h(0  )  x2t 


Solving  for  x.,  k  =  1,  2,  3 


n-1 
*1 


1  +  ki,Q  )t 


(120) 


.   k(en)x"-1  t  / 


•  .  /nnN  n-1  , 

,    k(8  )  x,    t 
r.    n-1   ,  /„n\4.  <"  n_l       y   1 
x  =  x    +  h(Q  )t  (x    +  - 


(8  )  xl    t     1- 

\         )     /   l+h(8n)t 

1  +  k(en)t     L 


(121) 
(122) 


The  problem  is  to  maximize  x  with  the  starting  conditions, 


(123) 
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Applying  the  maximum  principle  algorithm  and  referring  to 
Figure  4,  the  following  set  of  equations  are  obtained  from  equations 
(27).  (28).  (36)  and  (57): 


n_x     *£      t  k(en)  (,»  +  h(en)  A) 


1    l  +  k(en)t   (l  +  k(en)t)   (l  +  h(en)t) 

n 

1 


zn-l  _  4  -  "^)*  J 
2      1  +  h(On)t 


r.-l    n 
z3   "  z3 


The  boundary  conditions  for  z.  are 

l 


M      [1.  i-3 


z.    • 


1         0,   i/3 


Thus 


"3 


1  (124) 

■5  *  k(On)t 


.a-1       _2_ 

1   +   h(9n)t 

n_!         [l   +  h(Qn)tj      z?  +  k(9n)t     [,°  +   h(Qn)tj 

1  [l  +  h(en)t]     Pi  +  k(en)tj 


z2    -  ~ — rrzr  (125) 


(126) 


and 


i  ♦  k(en)t     [  2  J       2  i  +  k(on)t  L  J 

(127) 


The  following  substitutions  are  made  for  simplification: 
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P(8»)  =   "('"It 
1  +  k(on)t 


1  -  F(un)  =  — L 


n,    _  h(On)t 


G(0  . 


1  +  h(9n)t 


1  +  k(6n)t 


1  -  G(Gn) 


Fi(0  .  i_Fl£l 


d  0 


1  +  h(0")t 
(en  _d^l 


d  G 


Hn  -  x^-1  .J  [  i  -  P(en)J  ♦  [  x-1  +  p(en)  x^Vjfz^i  -  G(en)J  +  tG(en)) 


n-1 


The  optimal  decision  6  is  at 


9n    :    ■&£  .   0 


(128) 


x^z^ce11)  ♦  f1(ob)^"1  t  [.».[i  -  G(sn)J  +  to(enj}+ 

{x^1  +  Pie^-1  t]  [tGl(en)  -  ^(q11)]  .  o  (129) 


This  problem  does  not  yield  the  analytic  solution  as  did  the 
first  two  problems.   The  general  working  procedure  stated  before  should 
be  applied. 

A  flow  chart  for  solving  this  problem  is  shown  in  Figure  13,  and 
a  Fortran  1620  comruter  program  for  the  three-stage  system  is  included 
in  Appendix  C.   It  was  found  to  be  more  convenient  to  use  the 
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Start 

l 

Read    in     T,X?    X°    X? 

i 

Compute  xP 

n=  1,2, NK 

k  =  l.2.3 

1 

Compute    z" 

A 

—  o-o  — o-  <jL 

i 

a  — en 

i 

■ 

Compute  Hn 

,i, 

H°  =~dt.n1 

NO 

A  Yes    . 

Hn-~  Jln 

No 

A 

en+-A-^-en 

- — 1 

A 

» 

en>  b  ? 

A  Yes 

n-t-  I— n 

i 

No 

n   =    M  7 

J/  Yes 

Output  ^n    corresponding 
Qn  ax? 

i 

No 

Is  successive  e"x£  close 
enough? 

A  Yes 

Stop 

Fig.  13.      Flow    chart      for     maximum     principle 
solution      of      biochemical      reaction 
problem. 


,2 


expression  H     m  Kax  than 


*r 


0  to  solve  the  optimum  8  in  computation 


with  the  use  of  a  digital  computer. 

With  the  given  values  of  the  initial  concentrations,  x  =  1, 
x„  -   0.3,  x  =  0,  a  set  of  G  ,  8  and  8  is  arbitrarily  chosen, 
x  is  then  calculated  by  equations  (120),  (121)  and  (122)  and  z  by 
equations  (125)  slld  ( 126} .   Substituting  these  values  into  equation 
(127),  a  new  set  of  8  ,  8  and  8  is  obtained,   This  new  set  of 

decisions  is  employed  to  calculate  the  x  and  z  and  again  a  set  of 

12      3 
8  ,  8  and  6  .   The  procedure  is  stopped  when  two  consecutive  sets 

of  decisions  are  sufficiently  close.   Table  3  shows  the  detail  of 

calculation. 

The  machine  calculation  from  the  first  guess  to  the  final  answer 

required  less  than  half  an  hour  with  8  ranging  from  1.0  to  k.5   with 

an  increment  0.1.   Table  3  shows  that  *>  and  x  converge  very  rapidly, 

and  the  results  of  two  different  first  guesses  are  iuentical.  The 

computer  used  was  an  IBM  Fortran  1620  with  20,000  memories. 

The  dynamic  programming  solution  for  tku   same  given  initial 

12         3  ^ 

concentrations  were  8  =  2.0,  8  =  3.2,  Qy   =  ^.0,  and  x^  =  O.323I 

which  agree  very  closely  with  the  maximum  principle  solution.  A 

detailed  comparison  of  the  two  approaches  will  be  presented  later. 


Denbigh's  Reaction  System 
Denbigh(26)  propesed  the  following  reaction  systen 


(A  +  B) 
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TABLE  ^ 
Calculation  according  to  the  maximum  principle  algorithm 
Initial  concentrations  x1  =  1,  x  =  0.3,  x  =  0 

.   First  guessing  6  =  2.5,  82  =  3.0,  Q3  =  3.5 


First  Calculation 
Second  Calculation 
Third  Calculation 
Fourth  Calculation 
Fifth  Calculation 
Sixth  Calculation 


1         2         3 

First  Guessing  Q  =  4.0,  0  .  4.0,  6^  >  4.0 

First  Calculation 
Second  Calculation 
Third  Calculation 
Fourth  Calculation 
Fifth  Calculation 
Sixth  Calculation 
Seventh  Calculation 


e1 

e2 
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x3 
3 

2.7 

3-3 

4.0 

O.32O6 

2.6 

3.1 

3-9 

0.3208 

2.7 

3.2 

4.0 

0.3210 

2.6 

3.1 

4.0 

0.3210 

2.6 

3.2 

4.0 

C.3211 

2.6 

3-2 

4.0 

0.3211 

2.5 

2.8 

3.5 

O.3126 

2.8 

3-3 

4.1 

O.3196 

2.6 

3.1 

3.9 

O.32O8 

2-7 

3.2 

4.0 

O.32IO 

2.6 

3.1 

4.0 

0.3210 

2.6 

3-2 

4.0 

O.3211 

2.6 

3.2 

4.0 

O.3211 

6>r 


and  demonstrated  that  a  great  improvement  in  yield  of  the  desired 
product  Y  can  result  when  optimal  policies  are  adopted  for  a  stage- 
wise  reactor  system.   X  is  an  intermediate  and  P  and  Q  are  the  pro- 
ducts of  side  reactions.   A  and  B  are  the  raw  material  where  B  is  cheap 
and  excess  amounts  are  available  for  reaction.   This  system  comprises 
a  general  class  of  reactions  which  are  encountered  in  the  manufacturing 
chemical  industry  (26).  The  reactions  are  taken  to  be  first  order 
regarding  the  concentration  of  A  and  X.   The  reaction  rate  constants 
of  the  four  reactions  are  of  the  usual  Arrehnius  form 
k.  =  a.  exp  /-E./rt)   ,        i  =  1,  2,  3,  k 

Four  cases,  depending  on  the  relative  magnitude  of  the  reaction 
rate3,  arise  (6) . 

1.  E.  \E  ,  E ,>E,  :   It  can  he  seen  that  a  conversion  of  nearly 
all  A  to  Y  is  expected  with  a  sufficiently  high  operating 
tenperature. 

2.  3.<  E  ,  E  <E,:   For  this  case,  the  lowest  possible 
temperature  will  give  the  best  yield. 

3.  E.>  E  ,  E  <  E.  :   In  the  early  stages  of  the  reaction  the 
temperature  should  be  high  to  promote  the  production  of  X. 
As  soon  as  the  quantity  of  X  becomes  considerable,  the 
temperature  should  be  lowered  to  prevent  the  side  reaction. 
Hence,  a  decreasing  sequence  of  temperature  is  desirable. 

k.     E^  <  E  ,  '  E  y  E,  :   In  this  case  the  same  argument  calls  for 
an  increasing  temperature  sequence. 
The  first  two  cases  are  not  of  much  interest.   The  third  and  fourth 
cases  are  considered. 
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Let  w,  x  and  y  denote  the  concentrations  of  A,  X  and  Y  respec- 
tively and  'the  residence  time  in  each  stage.  The  material  'balance 
around  the  n-th  stage,  referring  to  Figure  2,  yields 


w  .  =  w  )   1  +  C      (k.  +  k„  ) 

n+1    n  (.     n  v  In  2n' 

x  .  =  -  w  P     k,   +  x  il  +  f      (k,   +  k,  ) 

n+1      n  >  n  In    n  ]  n  v  3n    <fn' 


yn+l  "  *  xn  e  n  KJn  +  yn 
Using  the  values  of  the  specific  reaction  rates  given  by  Denbigh  (26) 


—  =  10  exp  (-  -^-j-) 


kl 


*    in"*      /3000, 
—  .  10   exp  (^-) 


and  making  the  following  substitution 
k„ 


8  -  ^  .  10*  exp  (-  ^f  ) 


(130) 


the  material  balance  around  the  n-th  stage  of  a  system  consisting  of 
K  completely  stirred  tank  reactors  becomes 

(131) 

(132) 
(133) 


w   ,  =  fl  t  t  (1  +  8  )  I  w 
n+1    (.     n^     n'  \     n 


x  ,  -  -t  w  +   1  +  0.01  t   9"   (1  +  9  )  \   x 
n+1     n  n    )         n  n       n  \  n 


'n+1 


0.01  t  x  +  y 

n  n   ''n 
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Since  these  equations  are  homogeneous  in  the  concentrations,  the 
ratio  oi  the  state  variables  may  be  employed  (6)  to  reduce  the 
dimensionality  of  the  problem.  Aris  (6)  solved  this  problem  by 
means  of  the  cynamic  programming  algorithm,  proceeding  in  the  same 
way  as  done  in  solving  the  biochemical  reaction  problem.  When  he 
encountered  the  difficulty  of  the  individual  state  variables  being 
too  largo,  the  reciprocals  of  the  state  variables  were  employed  to 
improve  the  accuracy  of  interpolation  (6).  His  results  are  listed 
in  Table  *f. 

The  remainder  of  this  section  describes  the  solution  of  the 
problem  by  use  of  the  maximum  principle  algorithm.   The  results  from 
both  methods  are  also  compared. 

Adopting  the  notations  of  the  maximum  principle  algorithm 
and  referring  to  Figure  k,    equations  (151),  (152)  and  (155)  become 


x1   .  x  "J_  1  +  t  (I  +  0  )J 

n-1     n  n   f,    A  .  ,n,.        .n>  /.n]  n 
2   =      1  +  )   *      t  (1  +  6  )/Q  >  x„ 


n-1  n  n    n 

e,   ■  -  o.oi  t  x,  +  x. 
5  2    5 


Solving  for  x  ,  k  =  1,  2,  5, 


n-1 
**"  =  k n7  t1^ 

i  +  t  (i  +  e  ) 


n-1 


+  t 


1 


1  +  t"(l  +  a") 


i  +  o.oi  tn(i  +  Qn)/en 
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(U5) 


n    n-1 


n-1    ,  n 
K2   +  t 


3    3 


x'    +  0.01  t 


n-1 


1  +  t"(l  +  Qn) 


i  +  o.oi  tn(i:  h  on)  /e 


(136) 


The  problem  is  to  maximize  x  with  the  starting  conditions 
X°  =  0 


x°.o 


(137) 


Applying  the  maximum  principle  algorithm  the  following  set  of 
equations  corresponding  to  equations  (27),  (28),  (36)  and  (37)  is 
obtained 


n-1 


1 


z2 


1        "    1   +   tn(l    t   0")  [l   +   t"(l   +   Qn)J       [l    +    0.01    t"(l    +  -i)] 

0.01    z?(tn)2 
+  2_ . 

[1  +  tn(i  +  e11)]     h  +  0.01  tn(i  +  \)1 


n-1 


1   +   0 


,01  tn(l   +  ±)        1  +  0.01  tn(l  +  i) 


n-1  n 

z3       "  z3 


The  boundary  conditions  for  2.  are 


N    I  1.  i  -  3 
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z. 

'  0,  i  .  1,2 

Thus 

z3  "  X  n  -  1.  2 »  (138) 

n 

-°-1  .     8 ,    o.oi  tn ,   . 

2    i  +  o.oi  tn(i  +  -L)   i  +  o.oi  tn(i  +  X  )      U39; 
e  en 


z"  +  " 

n-l       zl  *  z 

zl 


n  _n 
2 


i  +  tn(i  ♦  Qn)   [i  +  tn(i  +  gn)J  [i  +  o.oi  tn(i  +  i)l 


_  2 

0.01  (tn) 


+  [i  I  tn(i  +  e11)]    [i  ^  o.oi  tn(i  I  ±)~1 


(HfO) 


and 


n-l  n 


xn_1  tn 
n-l      xl   t 


1  +  0.01  tn(l  +  -i-  )      <■  J     3 


i  +  t"(i  +  en)   /  i  +  o.oi  tn(i  ♦  -i.  ) 


e"  :   Hn  =  Max 


(1^1) 


The  general  procedure  for  solving  this  problem  is  exactly  the -same 
as  the  one  used  in  solving  the  biochemical  problem  in  the  last  sections. 
In  the  results  presented  in  Table  5,  the  notations  are  changed  slightly 
to  match  those  in  Table  k   for  the  convenience  of  comparison.  The 
Fortran  1620  programs  for  two  and  three-stage  problems  are  included 


TABLE  k 
Aris's  results  lor  Denbigh's  reactions  (6) 


Number  of 
Reactors  Temperature       "Holding  Tines"      Yield 

W  ?3    *a    Ix      tj   t      t]_         (>) 


With  no  restrictions 

1  326 


25 


2               280   =>o      Zf.o  o°  57. 4 

J  270    320   00      1.6   0.8  tx>  66.3 

With  restrictions  on  temperature  and  holding  time 

1              318     3^-0  2221 

2              276   ^A            2:82  2100  1+5.I 

3  260    288    39^  l.k      1.3  2100  k$.5 
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in  Appendices  J  and  E.  The  programs  can  be  easily  converted  into  a 
general  program  for  an  N  stage  system.  The  computing  time,  dependent 
linearly  on  the  number  of  stages,  is  almost  independent  of  the  number 
of  state  variables,   while  the  time  used  for  the  two-stage  system 
calculation  was  about  two  hours,  the  amount  of  time  would  largely  depend 
on  the  searching  range  and  searching  grid  of  the  two  variables.  The 
time  required  is  approximately  proportional  to  the  products  of  the 
number  of  increments  of  both  control  variables,  t   and  6  ,  used  in 
search.  Actually,  a  proper  choice  of  the  searching  grid  could  make 
the  calculation  of  a  three-stage  reaction  system  possible  in  half 
an  hour  without  the  sacrifice  of  too  much  accuracy. 

A  comparison  of  the  results  presented  in  Tables  k   and  5  shots 
that  they  are  very  close  in  the  case  without  constraints  on  the 
control  variables.   But  for  the  case  with  restrictions,  the  results 
do  show  sizable  differences  for  the  case  of  the  three-sta«;e  system. 
The  yields  depend  largely  on  the  control  variables  of  the  last  stage, 
and  the  difference  between  the  control  variables  of  the  last  sti.ge 
and  the  rest  of  the  stages  is  much  greater  in  the  first  case  than  is 
that  of  the  second  case.   This  difference  means  that  the  control  var- 
iables of  the  first  two  stages  have  more  influence  on  the  yields  in 
the  restricted  case  than  in  the  case  without  restriction.   Because  of 
the  inevitable  interpolation  in  the  dynamic  programming  solution,  a 
greater  error  would  result  for  the  case  with  restriction  than  for  the 
case  without  restricxion. 
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Another  fact  which  may  uphold  the  .statements  above,   is   that,   if  the  values 
of  the  control  variables  given  by  Aris  for  the  three-stage  process  with  re- 
striction are  used  to  calculate  the  yield  directly  from  equations   (131),   (132) 
and  (133),   the  answer  is  53.83$  instead  of  h9.5%  as  given  by  Aris.     It  appears 
that,  for  this  problem,   the  maximum  principle  approach  not  only  makes  the 
calculation  simpler  but  also  gives  more  accurate  results  than  dynamic  pro- 
gramming. 

Problem  of  Growth  and  Predation 

Some  interesting  models  of  bacterial  growth  and  competition  have  been 
proposed  recently.     In  this  section  two  problems  of  this  kind,  formulated  by 
Aris   (27)  using  the  dynamic  programming,  are  reformulated  by  using  the  maxi- 
mum principle  algorithm.     These  two  problems,   especially  the  second  one, 
exhibit  the  superiority  of  the  maximum  principle  approach  over  that  of  dynamic 
programming. 

The  dynamic  programming  solution  is  not  repeated.     A  short  statement 
of  the  problem,  however,  is  provided,  adopting  notations  in  Figure  h. 

First  Problem.     A  population  of  creatures  x    preys  on  a  population  x 

and  both  classes  are  undesirable   (27).     A  pesticide  is  effective  against 

n-1 
x  ,  but  not  against  x  .     If  x        is  the  population  of  x    at  the  beginning 

-i 

of  the  treatment  in  the  n-th  time  interval  and  x  "  the  corresponding 

population  of  x„  in  the  same  time  interval,  the  growth  and  predation  can  be 

represented  by 

n    n-1   ,  ,.n.  n-1 
x  -  x,   =  h(8  )  x 

11  1 
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Table  £ 
Denbigh's  system  calculated  by  the  maximum  principle  algorithm 


Number  of     Temperature  "Holding  Times"       Yield 

Reactors 

00        Tj   T2    Tx  tj    t2    tx  (/») 


With  no  restrictions 

1      326       «">  25 

2      279-9  o*       ^-55  c=o  57-99 

3  267.4  295  <x?  1.7  2.6  o"  66.55 

With  restrictions  on  temperature  and  holding  time 

1        318        340  22 

2        278.8  39^        3.85  2100  49.03 

3  265   292  39k  1.35  2.35  2100  5^.85 
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n         n-1         ,  ,ji,     n-1  ,  n.     n-1 

X2   "  X2       =       (      '  Xl  "^      ^  X2 


X*.-   [l  +  h(8n)]   x^1  (1^2) 

x^[l  +  m(en)j    xj"1     -k^x^1  (31,3) 


where  h,  k  and  m  are  functions  of  the  dosage  of  the  pesticide,  6,  which  is 
added  at  a  definite   time  interval-stage. 

The  problem  is  to  minimize  the  linear  combination  function  of  x    and  x 
at  the  end  of  the  pesticide  process,  with  the  initial  conditions  x     »  a, 
x     ■  b.     That  is  to  minimize   (ex    +  dx  ),  c  and  d  being  constants. 

Defining  a  new  state  variable 

n  n  j.  j  n 

x,  =  ex     +  dx„ 

3  12 


x^/o  [h  (8n)  +  i  j  -  dk(en)  I  +  d  [i  +  m(en)J  x^1 


The  problem  is  then  transformed  to  minimize  x  ,  with  the  starting  con- 

3 
ditions 

0 
x1-a 

x°-b 

0 
x  =  ac  +  bd 
3 

Applying  the  maximum  principle  algorithm  (5),  the  following  set  of 
equations  corresponding  to  equations  (27),  (28),  (36)  and  (37),  is  obtained: 
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"-1  =  fh(on)  *l}%l-  z\(en)  ♦  ■£  £  ( Me")  -  l) 
(e^  +  ^d  [i  +  m(en)J 


dk(8  ) 


1  +  m( 


n-1 
z    ■ 
2 

n-l   n 
z    =0 

3 


The  boundary  conditions  for  z_  are 


l,i  -  3 


10,1  -  1,2 


Thus 


z11'1  -  /h(en)  ♦  i 


n-l 


1  +  m(8n) 


n 


and 


Hn  -  ,J  £  ♦  h(e")  Jx""1  ♦  .jA.  *  m(9n)) x^1  -  k(en)  x^" 
H*  -  x^h^)  ♦  l}-  dk(8N)j  +  d  (l  ♦  m(8Ni)  x/-l 


(U*6) 

(11*7) 
(11*8) 

(11*9) 
(150) 
(1*1) 


8"  :  H"  -  Min 

n  -  1,2, ...,N 
The  general  procedure  for  solving  this  problem  is  to  follow  the  procedure 
stated  in  the  maximum  principle  section  by  using  equations  (11*2)  through 

(151). 

An  alternate  algorithm  (28,29)  to  solve  this  problem  without  introduc- 
ing the  third  state  variable  is  to  use  the  general  working  scheme  in  con- 
junction with  equations  (li*2),  (ll»3),  (11*7)  and  (1U8)  and  the  following 
boundary  conditions  instead  of  equations  (ll+5)  and  (ll*6): 
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0 
x  =  a 

1 

x„  =  b 


(152) 


N 

z  ■  c 
1 

z  «  d 
2 


Second  Problem.  A  nutrient,  x  ,  is  required  by  two  forms  of  an  organism, 
x  and  x  if  they  are  to  remain  viable,  x  is  not  capable  of  reproduction 
but  x  may  be  consumed  to  form  a  product  x  and  this  product  enhances  the 
reproductive  process.  If  the  control  action  is  taken  at  an  equal  time  interval, 
ssy  t,  the  following  relationship  among  x  can  be  deduced  by  a  simple  biolo- 
gical model: 

1.  the  spontaneous  change  from  x  to  x 

2    .1 

n-1    n     nt  ,  .  , 

X2   "  X2  "  ""V  (1?3) 

2.  the  formation  of  x    from  x  ,   the  reproduction  of  a  fraction,    (1  -  0n), 
of  itself  and  the  use  of  the  remaining  fraction  Sn  to  form  x 

x^-x^ax\+(l-en)x^t.eN^  (15U) 

3.  the  formation  of  x  and  its  use  to  affect  the  reproduction  con- 
stant x? 

k 
x^-x^be^t-^t  055, 

It,   the  effect  on  the  reproduction  constant  x11 

h 


n-1 
x,   -  xi 


?-«5t 


n    -  \  "  c5t  d56) 


76 


5.     the  use  of  the  nutrient  x    to  maintain  the  life  of  x    and  x  , 

xn-1  -  x"  =  -(dxn  +  lccn)t  (157) 

5  5  2  1 

where  a,  b,  c,  d  and  k  are  constants. 

From  these  relationships,  the  transformation  functions  can  be  derived 

from  equations  (153),  (155),  (156),  and  (157)  as  follows: 


n-1 


Vf\*)-4-   ^ 


T  (x^1,   0n)  -  x"  - 


at 

n-1  n  n 

x         -bOxt 

J 1_ 

1  -  ct 


1  (x11-1   en)  =  xn  =  xn-1  -  (    ct 


■  n-1 

1  .       tdx 

n-1      n.         n        n-1  2 

Vxk  ' e } "  x5  *  x5    +  rrn 


)   (x3       -  bO  x  t) 


tkx 


(158) 
(159) 
(160) 

(161) 


Combining  equations   (15U),    (158)  and  (160)     . 

5-S*^    (1  -  en)tenbt(xn)2  +  xn  h  -  At  t  (1  •-  en)t(xn"1  -   JSLy) 


n-1 


n-1 


„  ,       atx 
n-1  o 

-  x         + 


1  -  at 


(162) 


If  the  constants  are  available,   equation   (162),   which  is  a  quadratic 
equation,   can  be  solved.     Let  its  solution  be 


and  also  let 
D3 


n-1 


ctx 
n-1  3 


CrrSE>ZTibt»a*(i  -  e11)  +  1  -  e\t  ♦  (1  -  en)t(x^  -  r^)  (l63a) 
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Then  differentiating  equation  (162)  with  respect  to  x  ,  x , 

.  for  Zr\     ,    fl     ,  Zh  }h     ,     and 

Z.    n-l   -  n-l  n-l     n-l 

Jx          0  x  3  x  ^x 

1      2  3      k 


n-l       n-l 
X3     >\    ' 

and  x      ,  and  so] 

n-l 

?XJ 

T 
?    1 

-  D 

3     "-1 

<?   X 

1 

2 

atD 
1  -  at 

?\ 

ct(l  -  9n)tT  D 

n-l 
?x3 

1  -  ct 

(163) 


(100 


(165) 


i-^j  -  -  t1(i  -  en)tD  (166) 


7)     n-1 


(167) 


~n 


The  objective  is  to  choose  optimal  0  to  maximize  the  final  concentration 

N 
of  x  with  the  specified  initial  concentrations  of  the  x  and  with  the  exhaus- 
tion of  the  nutrient  as  the  end  point  of  the  process.  If  en  is  small,  the 

reproduction  constant  x  remains  small;  but,  if  it  is  large,  too  much  effort 
h 

is  put  into  enhancing  the  reproductive  process  and  too  little  x  is  left  to 
reproduce.  There  is  the  optical  decision  for  9n  to  maximize  the  production 
of  x  . 

The  problem  is  then  to  maximize  x  with  the  boundary  conditions 
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(168) 


and  the  constraint: 

0"  :  o^on^i  (169) 

This  is  a  fixed  end  point  problem.     Employing  equation  (27)  and  com- 
bining it  with  equations   (163)  through  (167) 

zl"    '  z5~z$  n  "  1»2,. ..,»  (170) 


ktz5   f  (171) 


n    n 
atDbte    z 

1  "  at     '   1  "  at     +   U  -  ct)(l  -at) 


atDcte  nbtz?  f  -} 

_  , i*  +  ZJ       dt  ktatD   /  (172} 

(i  -  at)(i  -  ct)     zs|  r^TE  -  rrsW  (  '  ' 

n-1         ot(l  -  On)tT  Dz"        f       Onbtct( 

z3    =  r^t^1  ♦  i T 


;(1  -  0")tT  D 
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ctz 


ktz^ 


1  -  ot       1  -  ct       1  -  ot 


/ 


z?-1  =  -Tn  (1  -  9n)tD     z?  -  J +  zM  +  __ 

h  1  /  1      i  _  ct       5 

with  the  boundary  condition  for  z. 


;tenbtz[| 


ct 


l,i  =  1 


zi  >,i/l,5 


1°' 


Thus  the  Hamiltonian  expression  is, 

Hn  =  z^  ♦  z^  +  z^  +  fa   ♦  zf$ 


(173) 


(171* ) 


(175) 


(176) 


end 


8"  :  H" 


Max 


(177) 


z^  is  treated  as  a  parameter,  like  the  A  in  the  sludge  problem. 
Assuming  a  value  of  z-,  the  problem  is  solved  by  the  general  working  pro- 
cedure without  recycle,  using  equations  (158)  through  (161)  and  (168) 
through  (177)  and  obtaining  the  optimal  sequence  of  decisions  Y  8/  for 

this  particular  z  .  The  value  of  z^  is  varied  until  the  calculated  xH 
5  5  5 

is  zero.  The  sequence  of  the  optimal  {  9n^  thus  obtained  is  the  answer 
for  the  problem. 

The  equation  derived  above  appear  to  be  lengthly  and  tedious,  but  the 
algorithm  is  straightforward  and  overcomes  the  difficulty  associated  with 
high  dimensionality,  encountered  in  the  dynamic  programming  solution,  and 
resulting  from  the  large  number  of  state  variables  involved  in  the  problem. 
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Problems  of  .Resource  Allocation 

In  this  short  section,   one  and  two-dimensional  allocation  processes 
formulated  by  Bellman  according  to  the  method  of  dynamic  programming   (U) 
are  reformulated  with  the  use  of  the  maximum  principle  algorithm.     A  gen- 
eral approach  to  the  problems  of  rnulti-dim6nsional  allocation  is  also 
suggested.     The  problem  of  allocation  of  nutrients  is  included  as  a  specific 
example. 

First,   the  case  of  allocating  one  type  of  resource  y  is    considered. 
The  resource  y  is  to  be  divided  into  N  independent  parts,  each  of  which  is 
called  an  activity.     There  is  a  utility  function  associated  with  each  activity. 
This  function  measures  the  dependence  of  the  return  from  this  activity  upon 
the  quantity  of  the  resource  allocated.     Let  h  be  the  utility  function  which 
is  the  net  profit  received  from  one  particular  activity.     If  the  allocation 
of  resource  to  the  n-th  activity  is  8  ,  the  objective  function  to  be  maxi- 
mized is   the  sum  of  the  interval  profits. 


=  Max 


r  n 


hn(en) 


with  the  constraints 


n-l 


8   >  0     ;         n-l,  2,...,N. 
The  maximum  principle  formulation  proceeds  as  follows: 

Defining 

n       .n,„n.  n-l 

x     =  h   (8   )  +  x 
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2    2 
with  the  boundary  conditions 

x°  =  0 

1 


4-0 


it  is  deduced  that 
N 


\=   2J    h  (8  ) 


n-1 

N 
N        <tt-       n 

x   -  2-    e   ■ r 

2        n-l 

N 
The  problem  is  then  transformed  into  that  of  maximizing  x  by  employing 

the  maximum  principle  algorithm  for  the  case  with  the  specified  end  point. 
Depending  on  the  nature  of  the  utility  function,  h(6),  the  solution  may 
be  analytical  or  semi-analytical  or  numerical. 

If  there  are  two  types  of  available  resources,  y  and  w,  to  be  allo- 
cated, the  objective  function  is 


0*=  Max  fh»(8»  e")| 


under  the  constraints 


21    -1 

n-l 


-  y 


Z      °2-» 


2 
n-l 
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8^270  n  =  1,2,..., N 

9n2;0  n  =  1.2,. ...N 

2 

Solution  by  the  maximum  principle  proceeds  in  the  manner  similar  to  the  case 
of  the  one-dimensional  allocations.     Defining 

n        n-1        n 
x     ■  x         +  0 
2         2  1 


n 

n- 

-1 

_n 

X 

X 

+ 

B„ 

3 

3 

2 

with  the  boundary  conditions 
x°  =  0  A 

0   . 
x  ■  0 
2 

x°  =  0 
3 

N 

x2-y 

N 
x  ■  w 

3 

N 
Now  the  problem  is  to  maximize  x  employing  the  maximum  principle  algorithm 

for  the  case  with  the  specified  end  point.  It  can  be  seen  that  an  n-dimen- 

sional  allocation  problem  with  n+1  state  variables  and  n  control  variables 

can  be  solved  similarily. 

If  the  function,  h(8)  is  quite  complicated,  the  method  of  dynamic 

programming  may  be  preferred  for  some  cases  of  one  and  two-dimensional  problems. 

However,  for  most  of  the  cases,  especially  when  analytical  solutions  are 
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permissible,  the  answers  are  much  easier  to  obtain  by  using  the  maximum  prin- 
ciple than  by  the  dynamic  programming  algorithm.  For  multi-dimensional  re- 
source allocation  problems,  the  dynamic  programming  approach  will  be  at  a 
impasse  because  of  the  exponential  increase  in  computational  labor;  the  maximum 
principle  method  may  be  the  only  choice. 

Suppose  that  there  are  of  N  species  of  food;  each  unit  of  the  individual 
food  species  contains  cf  units  of  calories  and  P  units  of  a  particular  nutrient 
content  which  is  to  be  maximized  in  the  diet.  Let  superscript  n  denote  the 
n-th  species  of  food  and  0n  be  the  units  of  each  food  species  to  be  allocated 
to  make  up  a  diet.  If  the  diet  requires  the  maximum  content  of  the  particular 
nutrient,  as  mentioned  before,  but  at  the  same  time,  requires  that  the  total 
units  of  calories  be  maintained  at  a  definite  value,  y,  the  objective  junction 
is, 


with  the  constraints 

H 

JT      ctnsn   -  y 
n-1 

9n->l     ,       n-1,   2,...,N. 
The  formulation  according  to  the  maximum  principle  proceeds  as  follows 
Defining 

^-J»uf*^  .  (i78) 

4  -  *rX  Hn  en  (179) 

with  the  boundary  conditions 
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(160) 


The  problem  is  transformed  into  that  of  maximizing  x  .     From  equation   (28), 


the  Hamiltonian  expression  is 


uH  _  „n,  n-1 

H     »  z   (x        +  Ji      8  )  +  z2(x  *<A      8  ) 


From  equation   (30) 


n         n-1 
z     •  z 
1         1 


(181) 


n        n-1 
z     «  z 
2         2 


«T  ■  z"  *  (182) 


Applying  the  boundary  conditions  of  z 


zN-l 


thus 


1 


■x  "  1    ,      n-1,  2,...,N  (183) 


^-^V^IW1'*/)''  (18U) 

8n  ;  H11  -  Max  (185) 


where 


e"2i 
Taking  z2  as  the  parameter,  the  solution  can  be  carried  out  by  the  gen- 
eral working  scheme  by  using  equations  (178),  (179),  (180),  (1810  and  (185). 
Zj  is  varied  until  the  boundary  conditions  are  met. 
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GENEliAL  DISCUSSION 

The  maximum  principle  and  the  dynamic  programming  are  generally  re- 
garded as  two  alternative  ways  of  solving  optimization  problems  of  non- 
linear nature  (30), 

The  principal  prerequisite  for  optimizing  a  process  by  means  of  the 
dynamic  programming  and  the  maximum  principle  algorithms  is  a  mathematical 
description  of  the  process  itself,  which  often  lags  far  behind  the  theory  of 
optimization.  The  principle  of  optimality,  however,  can  be  applied  to  any 
optimization  problems  of  stagewise  nature  irrespective  of  the  availability 
of  the  exact  mathematical  description  of  the  process.  The  directed  net 
work  problem  is  one  of  the  examples  for  which  the  transform  functions  are 
not  in  exact  functional  form. 

The  sludge  treatment  problem  is  an  example  for  which  the  analytical 
solution  can  be  obtained  either  by  the  maximum  principle  algorithm  or  by 
the  dynamic  programming  algorithm.  However,  for  many  of  the  problems, 
these  methods  are  not  always  practicable.  A  considerable  number  of  problems 
fall  into  a  class  for  which  semi -analytical  solutions  can  be  derived  by  use 
of  the  maximum  principle  algorithm,  as  in  the  case  of  the  s  tep  rocket  problem, 
but  it  is  not  so  when  the  dynamic  programming  algorithm  is  employed.  This 
resulting  in  semi -analytical  solution  is  one  of  the  advantages  in  using  the 
maximum  principle  algorithm.  The  semi-analytical  solution  reduces  computational 
labor  considerably.  For  some  problems  even  a  simple  hand  calculation  enables 
one  to  obtain  a  fairly  accurate  answer. 

The  majority  of  problems  belong  to  a  class  of  problems  for  which  no 
analytical  solutions  are  obtainable,  and  for  which  the  general  computing 
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scheme  must  be  employed.  The  last  three  problems  in  the  application  section 
belong  to  such  a  class  of  problems. 

Despite  outward  difference  the  basic  notions  behind  dynamic  programming 
and  the  maximum  principle  are  similar.  Under  most  circumstances  their  nu- 
merical answers  are  very  close  and  can  be  used  jointly.  Their  essential 
difference  is  a  tactical  one  (30).  While  dynamic  programming  starts  the 
investigation  by  searching  the  entire  grid  of  the  s  variables  at  one  stage, 
stores  this  grid  of  values,  and  proceeds  stage  by  stage,  the  maximum  prin- 
ciple starts  the  investigation  by  computing  one  optimum  path  along  the  n 
stages  and  then  proceeding  to  improve  this  optimum  path  based  on  the  values 
obtained  from  the  preceding  computation.  This  scheme  enables  the  maximum 
principle  to  deal  with  processes  in  which  the  optimum  conditions  at  any 
stage  can  be  disturbed  by  conditions  at  a  following  stage. 

The  dynamic  programming  algorithm  leads  immediately  to  a  neat  com- 
puter program  for  machine  calculation.  However,  since  it  must  store  the 
results  of  the  entire  grid  of  s  variables  at  each  stage,  the  computing 
time  may  become  considerable  and  require  a  large  capacity  of  computer 
memories.  Suppose  that  there  are  t  decisions  to  be  made  at  each  stage  for 
a  process  with  s  state  variables  at  each  stage  and  n  total  stages,  and 
also  suppose  that  a  grid  of  ten  for  each  variable — state  and  control — 
would  yield  the  results  with  sufficient  accuracy.  In  the  dynamic  program- 
ming formulation,  10*  calculations  would  be  required  at  each  stage  for  each 

s+t 
state  variable,  and  therefore,  a  total  of  10    calculations  would  be 

required  at  each  stage,  resulting  in  the  n  10    calculations  altogether. 

The  maximum  principle  formulation  would  need  n  10  calculations,  almost 

independent  of  the  increase  in  the  number  of  state  variables  s.  If  the 
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s  becomes  large,  dynamic  programming  may  require  much  more  calculation  than 
the  maximum  principle,  but  the  most  serious  difficulty  associated  with  the 
use  of  dynamic  programming  lies  in  the  exponential  increase  of  storage  space 
with  respect  to  s.  This  difficulty  is  one  of  the  reasons  that,  in  dynamic 
programming  formulation,  an  effort  is  always  made  to  reduce  the  number  of 
state  variables,  either  by  using  the  ratio  of  the  state  variables,  either 
by  using  the  ratio  of  the  state  variables  as  a  new  variable  or  by  using 
other  techniques.  Since  in  employing  dynamic  programming  the  computer 
program  must  be  written  by  double  indexing  the  variables,  memory  space 
is  easily  exhausted  by  an  increase  in  the  number  of  state  variables.  In 
programming  the  maximum  principle  algorithm  for  computation  only  one  index 
is  used  so  that  total  memory  space  increases  linearly  with  respect  to  s. 

The  error  introduced  due  to  the  interpolation  technique  employed  in 
the  dynamic  programming  algorithm  is  difficult  to  estimate.  If  the  grid  is 
sufficiently  small  and  the  function  has  no  flat  portion,  the  interpolation 
may  not  create  too  much  of  a  problem.  It  is  to  be  noted,  however,  that  an 
increase  in  grid  number  will  result  in  tremendous  increase  in  computing 
time  and  computer  memories.  Usually  a  balance  should  be  sought  in  light 
of  the  nature  of  the  transformation  functions. 

From  the  previous  statements  the  following  conclusion  may  be  drawn: 
If  an  appreciable  error  results  from  the  interpolation  and  the  number  of 
state  variables  is  large,  the  maximum  principle  approach  is  preferable.  If 
the  results  for  a  large  number  of  initial  values  of  the  state  variables 
are  desired,  it  would  be  too  tedious  to  run  the  maximum  principle  calcula- 
tions repeatedly.  In  such  instances  dynamic  programming  may  be  the  better 
choice.  If  the  results  for  only  a  few  initial  values  of  the  state  variables 
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are  wanted,  the  maximum  principle  approach  would  give  more  accurate  answers 
in  much  shorter  computing  time  than  dynamic  programming.  Generally  speak- 
ing, if  the  transformation  functions  are  complicated,  it  may  be  advisable 
to  work  with  dynamic  programming,  since  the  maximum  principle  algorithm 
requires  handling  of  equations  of  very  complex  form  which  may  easily  lead 
to  manipulation  error. 
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KCKENCLATURE 

a  =  a  given  constant 

a  =  a  given  constant 

a.  =  frequency  factor  in  Arrehnius  rate  equation 

a  =  -B  Mn  -  dL 
n  n       n 

A  =  chemical  specie 

b  =  a  given  constant 

b  =  a  given  constant 

b  =  1  -  c< 
n  n 

B  =  chemical  specie 

c  =  a  given  constant 

c  =  propellant  exhaust  velocity  of  stage  n 

C  =  chemical  specie 

d  =  a  given  constant 

,1  . .  _. 

a  =  a  given  constant 


d         =  w 
n  xn 


?(8  ; 


D  .    =  defined  by  equation  (163a) 

e  =  x/y 

E  =  sludge  alkalinity 

E.  =  the  activation  energy  of  the  chemical  reaction 

f  >  optimal  net  profit  of  the  time  interval  from  the 

present  stage  to  the  end  of  the  process 

f  (V)  =  minimum  weight  of  rocket  n  achieveing  velocity  V 


r-(on)t 
l  +  K(en)t 
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(8.}       .  <M^ 


f/x 


•O^v 


1  +  h(On)t 

1  d9n 

h         =  Optimal  net  profit  taking  A  as  a  parameter;  function 

of  control  variable;  chemical  reaction  rate  constant 
H         =  the  Hamiltonian 
i         =  node  number 
j         •  node  number;  an  indication  of  the  path  the  particle 

should  follow 
k        a  function  of  control  variable;  chemical  reaction 

constant;  a  given  constant 
m         =  function  of  the  control  variable,  total  grid  number 
K         =  ratio  of  initial  thrust  to  weight  of  rocket  n 
n         a  number  of  stages 
H         a  total  number  of  stages;  total  number  of  years  for 

a  resource  to  be  available;  lest  stage 
0         =  objective  function 
P         =  interal  profit;  chemical  specie 
q         =  flow  rate  of  the  feed  stream 
Q         =  chemical  specie 
GL-        =  total  amount  oi'  wash  water 
r         =  flow  rate  of  the  recycle  stream;  rate  of  reaction 

E         =  the  ideal  gas  law  constant;  flow  rate  of  wash  water 
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s         =  total  number  of  state  variables 

t         =  residence  time;  total  number  of  control  variable; 

control  variable 
T         =  transformation  function;  absolute  temperature 
u         =  x/y 
v         =  volumetric  flow  rate;  w/x 

V  =  reactor  volume;  design  velocity  of  rocket 

w         =  state  variable;  alkalinity  of  the  wash  water;  initial 

gross  weight  of  step  rocket 
w         =  that  portion  of  jettison  rocket  weight  which  is 

constant 

V  =  rocket  system  payload  weight 
L 

x  =  state  variable 

X  =   chemical  specie 

y  =  state  variable 

V  -  chemical  specie 

2         =  an  independent  variable;  a  variable  introduced  in  the 
basic  algorithm,  equation  (27) 

Greek  Letters 


^ 


=  a  given  constant 


P  =  a  given  constant 

C~\  =  a  given  constant 

X  =  Lagrange  multiplier;  relative  cost 

ft  =   Kronecker  delta,  defined  by  a  .      =  \     '    . 
a  '         "   ii   lo,  i 


©     =  residence  time;  control  variable 
[      =  residence  time 
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APPENDIX  A 


11 
12 


i4 


DIMENSi  ON  A(4)»3l4>»C(4)»D(4)»W<7v>4).VA<70»4) 

DIMENSICNWAI  7'-  .4)  >WB<  70.4)  »X(7C4)  >XA<70,4)  »XB(70>4)  .VI  70.4) 

F  C  R  K  A  T  (  E  1  0 .  4  ,  E  i  C  .  4  ,  E  1  C  •  4  ,  F.  1 0  .  4  ,  E 1  0  .  4  ,  E 1  0  .  4  ) 

R?ADliA(2)iA{9>,AU)»B{2)«BI3)t3l4) 

SEAD1  »C  t  2)  .C(3!iC(4)  .0(21.0(3)  .0(4) 

READ  1 ,? 

DC2N'=1  .70 

XA (M.2 I  =VN,*5C.O 

WA(M»2l *<  ! 1.C-AC2)  1  wpt-DI  2)  )  /  (EXP(-XA(M»2')  /C(2)  )-B(2)-A(2)  ) 

PUNCH I .WAIM.2 )  »XA(M  .2) 

D09I=3,4 

DC9N=1 .70 

KAiN.l  1=9999599999.9 

VN=l\i 

XIN.I )=3CC.C+VN*50.0 

V(N»I)=O.C 

X:.\.I-1)=X(N.I)-V(N»I) 

DC13K=1  ,7. 

XD=X(N. i-1 1-XAlM.l-l ) 

Ir  (XD--  .0.-)  11.14.12 

xo=-xj 

ir (XD~2d.C)14. 14.13 

L  J  )N  T  i  \  o  c 

.:  (N.I-1  l»X(N.I-l)*XA(M.I-l  1/XA1M.  1-1) 

•.'(.Nil  !  =  (  U.O-AI  1  >  )*W(N.I-1)+DU  )  )/<EXPl-V(,M.I  1 /C  (  I  )  ) -o ( I ) -A ( I i 

If (WAIN, I 1-WIN.I ) 16.6.5 

WA(N.I )=W(N. i  ) 

VA  ( N .  I  1  =V  ( H ,  I ) 

XA'(N»I)=X(M»I) 

WB<Nt  I-l,)=W(Ni  I -11 
Xo(N.I-i)=X(N« 1-1 1 
FA=X(N.l 1-50. C 
IF(  VI.N.i  )  -FA  17.9.9 

V1N.1 )=V(N.I 1+30.0 


PUNCH1.WAIN. I ) .VAiN.l I.XA(Ntl) »W8(N»I-1) >Xd( N.I-1) 
END 
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APPENDIX    S 


u  .  ij  >    <i.:\iE:TICS    PROBLEM 

D  X  i'icr>iS  I CN    'J  A  {  100  >4  )  »  0  t  100  »4 )  »r^l  i-ui'ti  »F  (  10 

KEAM.T 

DC    3    M= 1 » 1 CO 


Ust)  >  D  d  t  i  0 0  »  4  J 


DAi,'<;,  i  )  =VM*0.03 
r  A(.-i,l  )  -u 
DC6.\=2.A 

0^.       D       J  =  i  »  1  \J  U 

Vj=  J 

D  (  J  f  i\  J  -  V  J  *  w  •  u  J 

r  A(  J  »N)  =-U..VVV9c+iO 

P^i  »v> 

X=-P*"i  t9  .  i_-«-r-6  •  J 

D(  J»N-1  Js  (U(  JtN)*'ll.i;+Y*T)+Y*T  )  /(  i«0+X*Tl 
ci=X*T»(  D(  JiN  )  M  Y*T )/ (1.0+Y*1  !  )  /  ( 1.0+X*T  ! 
DC    12     «  =  1,1..'.! 
XD=D( Ji\-i)-OA<M»N-l ) 

i?  (  xn-.. . ;  1 1.,  ,ia  ,ii 

aC=-XJ 

IF { XD— ^ • o  i  3 )  i^>ib»12 

r  l«m\-.  /-ul  v  »|N— 1)*FA  I  ;•','.  i\-i  i  /  J  A  I  Mi,«-l  ) 
r  !  J  >  N )  =  ot  I-  (  J  > i\-l )  /  ( 1  .0  +  Y*l  J 
lr  lr  I  J.M-.FAI  J.N)  )    9  >9>b 
r  "  C  wi »  "n;  j  =  F  [  J  >  ,\i } 

*   IJ.N-l ) =D(J»N-1 1 
DA( j.N)=D(J.N) 
'l  F  [P-b.C  )2vt2C:  ,6 


3  IC'P 


i  U  >  i\  )  .  P  A  >  r  A  i  J  ,  N  )  i  Da  IJ  i  fv- 1  ) 
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■  y    -. . 


if? 


•j  i  ■'■'■:  hb  l  ON    6X  I  4J  >r  X  i  A-  )  iol»]  >  r  ( <»  ]  lOJlt)  ifulll  >  A  I  4  i  »  Y  !  4  J  i  U  ( <f  ) 
J  i  i-iiNS  i  JiN    Za  i  4  )  >ZY  I  *  )  »PA  (  4  )  iKa  (  4  )  >  f  i  4  )  >  KC  I  4  ) 

i    HvKMAT    ( '*.l.'«t»Ell>»<riEHJ«'t» EiO«4»fc'l^i^fi£10«'nE10i'») 

R  h  A  L<    I  >  T  »  P  (  2  1  »  P  {  jj  J  >  p  (  4  ) 

rtLnw      HAiii)hll),Ull)tJrt,J 

2    0  J    4    i\  -  2  >  4 

uX  (  N  )  =  -P  i  >N  1**2  +  9.0^  I  M)-b.O 

rXIN)=-3.  C*Pd\i  **2+10.0*P  (N)xii.O 

•j(N  )=1*oX  (N)'/(  1«0+T*0X  (N)  J 

F!.\1=T*FX(N>/(1  . 0+T*FX(N)  ) 

XI  N'  )  =  (  l.j-F(K)  )*X(N-1) 

Y  (  N  )  =  (  I  .  v-6I  N  )  I  *  (  Y  (  N-l  )  +F  (  N  !  *X  (  ft-  1  )  ) 

W  ( :■:  )  =\:  ( N-l ]+6l N)*<  Y(  N-l  )+F  (N)  *X  (  N-l  )  ) 
4    r  ?,  i  N  T     1  .  X  (  N  )  .  Y  t  M  )  .  W  (  i\  )  ,  P  (  N  )  >  u  X  ( t  i  )  >  r  X  ( i\  ) 

Z  X  (  U )  ■  -  .  ij 

ZY(4)=-.0 

ZXl  J  )=r  14  l*o(  4  ) 

£  Y  i  J  )  =  u  I  •+  ) 

4AI  i.  )-\  i.O-r  1  j  )  )  *iX(  3)i-r  i  J  )  "Itili)+Z.Y(3  !*(  i.0-o(  J  i  )  ) 

4  Y  t  Z  )  -  -  I  j  l  t  i  a  •  . :>  ~*o  t  j  )  i  '■'*-  Z  Y  v  3  ) 

Li  ..      O     \-  i  •  4 

:-'c  (  n  )  =-  vvs  9  v  y  ?9  s  9 .  99 

-    >':X  i  i\)  =-r-A  (  N)  »*2  +  9tO*P'A  (N1-6.C 

rXi  it  >"-'■>:  ^*PA(  N)**2  +  10.0*PA(  N.)  1-12.0 
rlN)-=i*rXlM/(i.O+T*hX(N)l 

o  i  N  )  =  1  * uX  (  N  )  /  I  i  .  i)+T*i>X  (  N  )  ) 

rJIM  =  i  ■■(  iu.O— 6.0*PA(N)  )  /  i  1.0  +  T*TXIN)  )**2 

ULilN)si*I>.j-iiCIVAlNI  )/( l.uf  T*(jXli\J  )  **z 

it  u  i  N  )  =A  i  .N  -1 )  *iX  (  N  )  *  (  1  .  0-r  (  N  )  )  tW  (  N-  i  I 

.<  =  i  N  i  =  kb  (  .N  )+  1  Y  I  i\-i  )  +r  (M*X  I  N-l  )  )  *  I  blNI+ZY  i  M*  I  i.U-u  I  i*  j  .  ) 

Ir  (,-.::  h\i-KClN)  )  10.10.9 
■<     iv  .'..  (  l\  /  =  :\  -j  t  M  ) 

i  f     .- .-,  (  k  )=r>A  (  i\  )-rPAP  • 
i  F  (  PA  (  N')-4.5  )  3  .  6  .6 
'>    P.JwCn    1  tSC(N)  tPINJ 

bL    TO    i 
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•  u+£  i  i\  )*  (  1 «  0+P  (  N  }  )  ) 
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A  brief  introduction  to  the  principle  of  optimality  is  given 
at  the  beginning  by  an  illustration  of  it  using  the  directed  net 
work  problem.   Brief  summaries  of  dynamic  programming  and  the 
maximum  principle  algorithms  follow. 

In  the  application  section,  six  problems  are  solved  and  pre- 
sented in  order  oi  complexity.   The  sludge  treatment  problem  is  a 
fairly  simple  case  for  which  identical  analytical  solutions  could 
be  obtained  by  both  dynamic  programming  and  the  maximum  principle 
algorithms.  The  step  rocket  problem  resulted  in  a  semi-analytical 
solution  by  the  use  of  the  end  specification  algorithm  of  the  maximum 
principle.  Regarding  computational,  this  method  is  often  much  simpler 
than  dynamic  programming  solution.   The  third  problem  is  the  optimization 
of  a  general  consecutive  first  order  biocheir.ical  reaction  system  in  a 
series  oi  ideal  baokmix  reactors.   The  general  working  schemes  were 
employed.   Tie  numerical  solutions  obtained  by  both  methods  were 
compared  in  detail.   Denbigh's  reaction  system  was  solved  by  means  of 
the  maximum  principle  ai.d  the  answers  were  compared  with  those 
obtained  by  Aris  with  the  use  of  the  dynamic  programming  approach. 

The  section  on  application  also  includes  two  problems  of  growth 
and  predation.   The  solutions  of  these  two  problems,  especially  the 
second  one,  clearly  show  the  advantage  of  employing  the  maximum 
principle  approach  over  employing  the  dynamic  programming  method 
when  the  number  of  state  variables  is  large. 

The  last  section  gives  a  general  discussion  on  the  relative 
merits  of.  the  two  approaches. 


